Nektar++
Loading...
Searching...
No Matches
Expansion3D.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: Expansion3D.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 Expansion3D routines
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <LibUtilities/Foundations/Interp.h>
36#include <LibUtilities/Foundations/InterpCoeff.h>
37#include <LocalRegions/Expansion2D.h>
38#include <LocalRegions/Expansion3D.h>
39#include <LocalRegions/MatrixKey.h>
40#include <LocalRegions/QuadExp.h>
41#include <LocalRegions/TriExp.h>
42#include <SpatialDomains/Geometry3D.h>
43
44using namespace std;
45
46namespace Nektar::LocalRegions
47{
48// evaluate additional terms in HDG face. Note that this assumes that
49// edges are unpacked into local cartesian order.
50void Expansion3D::AddHDGHelmholtzFaceTerms(
51 const NekDouble tau, const int face, Array<OneD, NekDouble> &facePhys,
52 const StdRegions::VarCoeffMap &varcoeffs, Array<OneD, NekDouble> &outarray)
53{
54 ExpansionSharedPtr FaceExp = GetTraceExp(face);
55 int i, j, n;
56 int nquad_f = FaceExp->GetNumPoints(0) * FaceExp->GetNumPoints(1);
57 int order_f = FaceExp->GetNcoeffs();
58 int coordim = GetCoordim();
59 int ncoeffs = GetNcoeffs();
60 bool mmf = (varcoeffs.find(StdRegions::eVarCoeffMF1x) != varcoeffs.end());
61
62 Array<OneD, NekDouble> inval(nquad_f);
63 Array<OneD, NekDouble> outcoeff(order_f);
64 Array<OneD, NekDouble> tmpcoeff(ncoeffs);
65
66 const Array<OneD, const Array<OneD, NekDouble>> &normals =
67 GetTraceNormal(face);
68
69 DNekScalMat &invMass = *GetLocMatrix(StdRegions::eInvMass);
70
71 DNekVec Coeffs(ncoeffs, outarray, eWrapper);
72 DNekVec Tmpcoeff(ncoeffs, tmpcoeff, eWrapper);
73
74 IndexMapKey ikey(eFaceToElement, DetShapeType(), GetBasisNumModes(0),
75 GetBasisNumModes(1), GetBasisNumModes(2), face,
76 GetTraceOrient(face));
77 IndexMapValuesSharedPtr map = GetIndexMap(ikey);
78
79 StdRegions::MatrixType DerivType[3] = {StdRegions::eWeakDeriv0,
80 StdRegions::eWeakDeriv1,
81 StdRegions::eWeakDeriv2};
82
83 // @TODO Variable coefficients
84 /*
85 StdRegions::VarCoeffType VarCoeff[3] = {StdRegions::eVarCoeffD00,
86 StdRegions::eVarCoeffD11,
87 StdRegions::eVarCoeffD22};
88 Array<OneD, NekDouble> varcoeff_work(nquad_f);
89 StdRegions::VarCoeffMap::const_iterator x;
90 ///// @TODO: What direction to use here??
91 if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
92 {
93 GetPhysFaceVarCoeffsFromElement(face,FaceExp,x->second,varcoeff_work);
94 Vmath::Vmul(nquad_f,varcoeff_work,1,FaceExp->GetPhys(),1,FaceExp->UpdatePhys(),1);
95 }
96 */
97
98 //================================================================
99 // Add F = \tau <phi_i,in_phys>
100 // Fill face and take inner product
101 FaceExp->IProductWRTBase(facePhys, outcoeff);
102
103 for (i = 0; i < order_f; ++i)
104 {
105 outarray[(*map)[i].index] += (*map)[i].sign * tau * outcoeff[i];
106 }
107 //================================================================
108
109 //===============================================================
110 // Add -\sum_i D_i^T M^{-1} G_i + E_i M^{-1} G_i =
111 // \sum_i D_i M^{-1} G_i term
112
113 // Three independent direction
114 for (n = 0; n < coordim; ++n)
115 {
116 if (mmf)
117 {
118 StdRegions::VarCoeffMap Weight;
119 Weight[StdRegions::eVarCoeffMass] = GetMFMag(n, varcoeffs);
120
121 MatrixKey invMasskey(StdRegions::eInvMass, DetShapeType(), *this,
122 StdRegions::NullConstFactorMap, Weight);
123
124 invMass = *GetLocMatrix(invMasskey);
125
126 Array<OneD, NekDouble> ncdotMF_f =
127 GetnFacecdotMF(n, face, FaceExp, normals, varcoeffs);
128
129 Vmath::Vmul(nquad_f, ncdotMF_f, 1, facePhys, 1, inval, 1);
130 }
131 else
132 {
133 Vmath::Vmul(nquad_f, normals[n], 1, facePhys, 1, inval, 1);
134 }
135
136 NekDouble scale = invMass.Scale();
137 const NekDouble *data = invMass.GetRawPtr();
138
139 // @TODO Multiply by variable coefficients
140 // @TODO: Document this (probably not needed)
141 /*
142 StdRegions::VarCoeffMap::const_iterator x;
143 if ((x = varcoeffs.find(VarCoeff[n])) != varcoeffs.end())
144 {
145 GetPhysEdgeVarCoeffsFromElement(edge,FaceExp,x->second,varcoeff_work);
146 Vmath::Vmul(nquad_f,varcoeff_work,1,FaceExp->GetPhys(),1,FaceExp->UpdatePhys(),1);
147 }
148 */
149
150 FaceExp->IProductWRTBase(inval, outcoeff);
151
152 // M^{-1} G
153 for (i = 0; i < ncoeffs; ++i)
154 {
155 tmpcoeff[i] = 0;
156 for (j = 0; j < order_f; ++j)
157 {
158 tmpcoeff[i] += scale * data[i + (*map)[j].index * ncoeffs] *
159 (*map)[j].sign * outcoeff[j];
160 }
161 }
162
163 if (mmf)
164 {
165 StdRegions::VarCoeffMap VarCoeffDirDeriv;
166 VarCoeffDirDeriv[StdRegions::eVarCoeffMF] =
167 GetMF(n, coordim, varcoeffs);
168 VarCoeffDirDeriv[StdRegions::eVarCoeffMFDiv] =
169 GetMFDiv(n, varcoeffs);
170
171 MatrixKey Dmatkey(StdRegions::eWeakDirectionalDeriv, DetShapeType(),
172 *this, StdRegions::NullConstFactorMap,
173 VarCoeffDirDeriv);
174
175 DNekScalMat &Dmat = *GetLocMatrix(Dmatkey);
176
177 Coeffs = Coeffs + Dmat * Tmpcoeff;
178 }
179 else
180 {
181 DNekScalMat &Dmat = *GetLocMatrix(DerivType[n]);
182 Coeffs = Coeffs + Dmat * Tmpcoeff;
183 }
184
185 /*
186 if(varcoeffs.find(VarCoeff[n]) != varcoeffs.end())
187 {
188 MatrixKey mkey(DerivType[n], DetExpansionType(), *this,
189 StdRegions::NullConstFactorMap, varcoeffs); DNekScalMat &Dmat =
190 *GetLocMatrix(mkey); Coeffs = Coeffs + Dmat*Tmpcoeff;
191 }
192
193 else
194 {
195 DNekScalMat &Dmat = *GetLocMatrix(DerivType[n]);
196 Coeffs = Coeffs + Dmat*Tmpcoeff;
197 }
198 */
199 }
200}
201
202void Expansion3D::GetPhysFaceVarCoeffsFromElement(
203 const int face, ExpansionSharedPtr &FaceExp,
204 const Array<OneD, const NekDouble> &varcoeff,
205 Array<OneD, NekDouble> &outarray)
206{
207 Array<OneD, NekDouble> tmp(GetNcoeffs());
208 Array<OneD, NekDouble> facetmp(FaceExp->GetNcoeffs());
209
210 // FwdTrans varcoeffs
211 FwdTrans(varcoeff, tmp);
212
213 // Map to edge
214 Array<OneD, unsigned int> emap;
215 Array<OneD, int> sign;
216
217 GetTraceToElementMap(face, emap, sign, GetTraceOrient(face));
218
219 for (unsigned int i = 0; i < FaceExp->GetNcoeffs(); ++i)
220 {
221 facetmp[i] = tmp[emap[i]];
222 }
223
224 // BwdTrans
225 FaceExp->BwdTrans(facetmp, outarray);
226}
227
228/**
229 * Computes the C matrix entries due to the presence of the identity
230 * matrix in Eqn. 32.
231 */
232void Expansion3D::AddNormTraceInt(const int dir,
233 Array<OneD, const NekDouble> &inarray,
234 Array<OneD, ExpansionSharedPtr> &FaceExp,
235 Array<OneD, NekDouble> &outarray,
236 const StdRegions::VarCoeffMap &varcoeffs)
237{
238 int i, f, cnt;
239 int order_f, nquad_f;
240 int nfaces = GetNtraces();
241
242 cnt = 0;
243 for (f = 0; f < nfaces; ++f)
244 {
245 order_f = FaceExp[f]->GetNcoeffs();
246 nquad_f = FaceExp[f]->GetNumPoints(0) * FaceExp[f]->GetNumPoints(1);
247
248 const Array<OneD, const Array<OneD, NekDouble>> &normals =
249 GetTraceNormal(f);
250 Array<OneD, NekDouble> faceCoeffs(order_f);
251 Array<OneD, NekDouble> facePhys(nquad_f);
252
253 for (i = 0; i < order_f; ++i)
254 {
255 faceCoeffs[i] = inarray[i + cnt];
256 }
257 cnt += order_f;
258
259 FaceExp[f]->BwdTrans(faceCoeffs, facePhys);
260
261 // Multiply by variable coefficient
262 /// @TODO: Document this
263 // StdRegions::VarCoeffType VarCoeff[3] =
264 // {StdRegions::eVarCoeffD00,
265 // StdRegions::eVarCoeffD11,
266 // StdRegions::eVarCoeffD22};
267 // StdRegions::VarCoeffMap::const_iterator x;
268 // Array<OneD, NekDouble> varcoeff_work(nquad_e);
269
270 // if ((x = varcoeffs.find(VarCoeff[dir])) !=
271 // varcoeffs.end())
272 // {
273 // GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
274 // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
275 // }
276 StdRegions::VarCoeffMap::const_iterator x;
277 if ((x = varcoeffs.find(StdRegions::eVarCoeffMF1x)) != varcoeffs.end())
278 {
279 Array<OneD, NekDouble> ncdotMF_f =
280 GetnFacecdotMF(dir, f, FaceExp[f], normals, varcoeffs);
281
282 Vmath::Vmul(nquad_f, ncdotMF_f, 1, facePhys, 1, facePhys, 1);
283 }
284 else
285 {
286 Vmath::Vmul(nquad_f, normals[dir], 1, facePhys, 1, facePhys, 1);
287 }
288
289 AddFaceBoundaryInt(f, FaceExp[f], facePhys, outarray, varcoeffs);
290 }
291}
292
293// shorter version of the above (coefficients are already set for faces)
294void Expansion3D::AddNormTraceInt(
295 const int dir, Array<OneD, ExpansionSharedPtr> &FaceExp,
296 Array<OneD, Array<OneD, NekDouble>> &faceCoeffs,
297 Array<OneD, NekDouble> &outarray)
298{
299 int f;
300 int nquad_f;
301 int nfaces = GetNtraces();
302
303 for (f = 0; f < nfaces; ++f)
304 {
305 nquad_f = FaceExp[f]->GetNumPoints(0) * FaceExp[f]->GetNumPoints(1);
306
307 const Array<OneD, const Array<OneD, NekDouble>> &normals =
308 GetTraceNormal(f);
309 Array<OneD, NekDouble> facePhys(nquad_f);
310
311 FaceExp[f]->BwdTrans(faceCoeffs[f], facePhys);
312
313 Vmath::Vmul(nquad_f, normals[dir], 1, facePhys, 1, facePhys, 1);
314
315 AddFaceBoundaryInt(f, FaceExp[f], facePhys, outarray);
316 }
317}
318
319/**
320 * For a given face add the \tilde{F}_1j contributions
321 */
322void Expansion3D::AddFaceBoundaryInt(
323 const int face, ExpansionSharedPtr &FaceExp,
324 Array<OneD, NekDouble> &facePhys, Array<OneD, NekDouble> &outarray,
325 [[maybe_unused]] const StdRegions::VarCoeffMap &varcoeffs)
326{
327 int i;
328 int order_f = FaceExp->GetNcoeffs();
329 Array<OneD, NekDouble> coeff(order_f);
330
331 IndexMapKey ikey(eFaceToElement, DetShapeType(), GetBasisNumModes(0),
332 GetBasisNumModes(1), GetBasisNumModes(2), face,
333 GetTraceOrient(face));
334 IndexMapValuesSharedPtr map = GetIndexMap(ikey);
335
336 FaceExp->IProductWRTBase(facePhys, coeff);
337
338 // add data to out array
339 for (i = 0; i < order_f; ++i)
340 {
341 outarray[(*map)[i].index] += (*map)[i].sign * coeff[i];
342 }
343}
344
345/**
346 * @brief Align face orientation with the geometry orientation.
347 */
348void Expansion3D::SetFaceToGeomOrientation(const int face,
349 Array<OneD, NekDouble> &inout)
350{
351 int j, k;
352 int nface = GetTraceNcoeffs(face);
353 Array<OneD, NekDouble> f_in(nface);
354 Vmath::Vcopy(nface, &inout[0], 1, &f_in[0], 1);
355
356 // retreiving face to element map for standard face orientation and
357 // for actual face orientation
358 IndexMapKey ikey1(eFaceToElement, DetShapeType(), GetBasisNumModes(0),
359 GetBasisNumModes(1), GetBasisNumModes(2), face,
360 StdRegions::eDir1FwdDir1_Dir2FwdDir2);
361 IndexMapValuesSharedPtr map1 = GetIndexMap(ikey1);
362 IndexMapKey ikey2(eFaceToElement, DetShapeType(), GetBasisNumModes(0),
363 GetBasisNumModes(1), GetBasisNumModes(2), face,
364 GetTraceOrient(face));
365 IndexMapValuesSharedPtr map2 = GetIndexMap(ikey2);
366
367 ASSERTL1((*map1).size() == (*map2).size(),
368 "There is an error with the GetTraceToElementMap");
369
370 for (j = 0; j < (*map1).size(); ++j)
371 {
372 // j = index in the standard orientation
373 for (k = 0; k < (*map2).size(); ++k)
374 {
375 // k = index in the actual orientation
376 if ((*map1)[j].index == (*map2)[k].index && k != j)
377 {
378 inout[k] = f_in[j];
379 // checking if sign is changing
380 if ((*map1)[j].sign != (*map2)[k].sign)
381 {
382 inout[k] *= -1.0;
383 }
384 break;
385 }
386 }
387 }
388}
389
390/**
391 * @brief Align trace orientation with the geometry orientation.
392 */
393void Expansion3D::SetTraceToGeomOrientation(Array<OneD, NekDouble> &inout)
394{
395 int i, cnt = 0;
396 int nfaces = GetNtraces();
397
398 Array<OneD, NekDouble> f_tmp;
399
400 for (i = 0; i < nfaces; ++i)
401 {
402 SetFaceToGeomOrientation(i, f_tmp = inout + cnt);
403 cnt += GetTraceNcoeffs(i);
404 }
405}
406
407DNekScalMatSharedPtr Expansion3D::CreateMatrix(const MatrixKey &mkey)
408{
409 DNekScalMatSharedPtr returnval;
410 LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
411
412 ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,
413 "Geometric information is not set up");
414
415 switch (mkey.GetMatrixType())
416 {
417 case StdRegions::eMass:
418 {
419 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
420 mkey.GetNVarCoeff())
421 {
422 NekDouble one = 1.0;
423 DNekMatSharedPtr mat = GenMatrix(mkey);
424 returnval =
425 MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
426 }
427 else
428 {
429 NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
430 DNekMatSharedPtr mat = GetStdMatrix(mkey);
431 returnval =
432 MemoryManager<DNekScalMat>::AllocateSharedPtr(jac, mat);
433 }
434 }
435 break;
436 case StdRegions::eInvMass:
437 {
438 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
439 {
440 NekDouble one = 1.0;
441 StdRegions::StdMatrixKey masskey(StdRegions::eMass,
442 DetShapeType(), *this);
443 DNekMatSharedPtr mat = GenMatrix(masskey);
444 mat->Invert();
445 returnval =
446 MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
447 }
448 else
449 {
450 NekDouble fac = 1.0 / (m_metricinfo->GetJac(ptsKeys))[0];
451 DNekMatSharedPtr mat = GetStdMatrix(mkey);
452 returnval =
453 MemoryManager<DNekScalMat>::AllocateSharedPtr(fac, mat);
454 }
455 }
456 break;
457 case StdRegions::eWeakDeriv0:
458 case StdRegions::eWeakDeriv1:
459 case StdRegions::eWeakDeriv2:
460 {
461 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
462 mkey.GetNVarCoeff())
463 {
464 NekDouble one = 1.0;
465 DNekMatSharedPtr mat = GenMatrix(mkey);
466
467 returnval =
468 MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
469 }
470 else
471 {
472 NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
473 Array<TwoD, const NekDouble> df =
474 m_metricinfo->GetDerivFactors(ptsKeys);
475 int dir = 0;
476 if (mkey.GetMatrixType() == StdRegions::eWeakDeriv0)
477 {
478 dir = 0;
479 }
480 if (mkey.GetMatrixType() == StdRegions::eWeakDeriv1)
481 {
482 dir = 1;
483 }
484 if (mkey.GetMatrixType() == StdRegions::eWeakDeriv2)
485 {
486 dir = 2;
487 }
488
489 MatrixKey deriv0key(StdRegions::eWeakDeriv0,
490 mkey.GetShapeType(), *this);
491 MatrixKey deriv1key(StdRegions::eWeakDeriv1,
492 mkey.GetShapeType(), *this);
493 MatrixKey deriv2key(StdRegions::eWeakDeriv2,
494 mkey.GetShapeType(), *this);
495
496 DNekMat &deriv0 = *GetStdMatrix(deriv0key);
497 DNekMat &deriv1 = *GetStdMatrix(deriv1key);
498 DNekMat &deriv2 = *GetStdMatrix(deriv2key);
499
500 int rows = deriv0.GetRows();
501 int cols = deriv1.GetColumns();
502
503 DNekMatSharedPtr WeakDeriv =
504 MemoryManager<DNekMat>::AllocateSharedPtr(rows, cols);
505 (*WeakDeriv) = df[3 * dir][0] * deriv0 +
506 df[3 * dir + 1][0] * deriv1 +
507 df[3 * dir + 2][0] * deriv2;
508
509 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(
510 jac, WeakDeriv);
511 }
512 }
513 break;
514 case StdRegions::eLaplacian:
515 {
516 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
517 (mkey.HasVarCoeff(StdRegions::eVarCoeffLaplacian)) ||
518 (mkey.HasVarCoeff(StdRegions::eVarCoeffD00)) ||
519 (mkey.HasVarCoeff(StdRegions::eVarCoeffD01)) ||
520 (mkey.HasVarCoeff(StdRegions::eVarCoeffD10)) ||
521 (mkey.HasVarCoeff(StdRegions::eVarCoeffD02)) ||
522 (mkey.HasVarCoeff(StdRegions::eVarCoeffD20)) ||
523 (mkey.HasVarCoeff(StdRegions::eVarCoeffD11)) ||
524 (mkey.HasVarCoeff(StdRegions::eVarCoeffD12)) ||
525 (mkey.HasVarCoeff(StdRegions::eVarCoeffD21)) ||
526 (mkey.HasVarCoeff(StdRegions::eVarCoeffD22)) ||
527 (mkey.ConstFactorExists(StdRegions::eFactorSVVCutoffRatio)))
528 {
529 NekDouble one = 1.0;
530 DNekMatSharedPtr mat = GenMatrix(mkey);
531
532 returnval =
533 MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
534 }
535 else
536 {
537 MatrixKey lap00key(StdRegions::eLaplacian00,
538 mkey.GetShapeType(), *this);
539 MatrixKey lap01key(StdRegions::eLaplacian01,
540 mkey.GetShapeType(), *this);
541 MatrixKey lap02key(StdRegions::eLaplacian02,
542 mkey.GetShapeType(), *this);
543 MatrixKey lap11key(StdRegions::eLaplacian11,
544 mkey.GetShapeType(), *this);
545 MatrixKey lap12key(StdRegions::eLaplacian12,
546 mkey.GetShapeType(), *this);
547 MatrixKey lap22key(StdRegions::eLaplacian22,
548 mkey.GetShapeType(), *this);
549
550 DNekMat &lap00 = *GetStdMatrix(lap00key);
551 DNekMat &lap01 = *GetStdMatrix(lap01key);
552 DNekMat &lap02 = *GetStdMatrix(lap02key);
553 DNekMat &lap11 = *GetStdMatrix(lap11key);
554 DNekMat &lap12 = *GetStdMatrix(lap12key);
555 DNekMat &lap22 = *GetStdMatrix(lap22key);
556
557 NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
558 Array<TwoD, const NekDouble> gmat =
559 m_metricinfo->GetGmat(ptsKeys);
560
561 int rows = lap00.GetRows();
562 int cols = lap00.GetColumns();
563
564 DNekMatSharedPtr lap =
565 MemoryManager<DNekMat>::AllocateSharedPtr(rows, cols);
566
567 (*lap) = gmat[0][0] * lap00 + gmat[4][0] * lap11 +
568 gmat[8][0] * lap22 +
569 gmat[3][0] * (lap01 + Transpose(lap01)) +
570 gmat[6][0] * (lap02 + Transpose(lap02)) +
571 gmat[7][0] * (lap12 + Transpose(lap12));
572
573 returnval =
574 MemoryManager<DNekScalMat>::AllocateSharedPtr(jac, lap);
575 }
576 }
577 break;
578 case StdRegions::eHelmholtz:
579 {
580 NekDouble factor = mkey.GetConstFactor(StdRegions::eFactorLambda);
581 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
582 DNekScalMat &MassMat = *GetLocMatrix(masskey);
583 MatrixKey lapkey(StdRegions::eLaplacian, mkey.GetShapeType(), *this,
584 mkey.GetConstFactors(), mkey.GetVarCoeffs());
585 DNekScalMat &LapMat = *GetLocMatrix(lapkey);
586
587 int rows = LapMat.GetRows();
588 int cols = LapMat.GetColumns();
589
590 DNekMatSharedPtr helm =
591 MemoryManager<DNekMat>::AllocateSharedPtr(rows, cols);
592
593 NekDouble one = 1.0;
594 (*helm) = LapMat + factor * MassMat;
595
596 returnval =
597 MemoryManager<DNekScalMat>::AllocateSharedPtr(one, helm);
598 }
599 break;
600 case StdRegions::eHelmholtzGJP:
601 {
602 MatrixKey helmkey(mkey, StdRegions::eHelmholtz);
603 DNekScalMat &HelmMat = *GetLocMatrix(helmkey);
604
605 // Generate a local copy of traceMat
606 MatrixKey key(mkey, StdRegions::eNormDerivOnTrace);
607 DNekMatSharedPtr NDTraceMat = Expansion3D::v_GenMatrix(key);
608
609 ASSERTL1(mkey.ConstFactorExists(StdRegions::eFactorGJP),
610 "Need to specify eFactorGJP to construct "
611 "a HelmholtzGJP matrix");
612
613 NekDouble factor = mkey.GetConstFactor(StdRegions::eFactorGJP);
614
615 factor /= HelmMat.Scale();
616
617 int ntot = HelmMat.GetRows() * HelmMat.GetColumns();
618
619 Vmath::Svtvp(ntot, factor, &NDTraceMat->GetPtr()[0], 1,
620 HelmMat.GetRawPtr(), 1, &NDTraceMat->GetPtr()[0], 1);
621
622 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(
623 HelmMat.Scale(), NDTraceMat);
624 }
625 break;
626 case StdRegions::eLinearAdvectionReaction:
627 {
628 NekDouble lambda = mkey.GetConstFactor(StdRegions::eFactorLambda);
629
630 // Construct mass matrix
631 // Check for mass-specific varcoeffs to avoid unncessary
632 // re-computation of the elemental matrix every time step
633 StdRegions::VarCoeffMap massVarcoeffs = StdRegions::NullVarCoeffMap;
634 if (mkey.HasVarCoeff(StdRegions::eVarCoeffMass))
635 {
636 massVarcoeffs[StdRegions::eVarCoeffMass] =
637 mkey.GetVarCoeff(StdRegions::eVarCoeffMass);
638 }
639 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this,
640 mkey.GetConstFactors(), massVarcoeffs);
641 DNekScalMat &MassMat = *GetLocMatrix(masskey);
642
643 // Construct advection matrix
644 // Check for varcoeffs not required;
645 // assume advection velocity is always time-dependent
646 MatrixKey advkey(mkey, StdRegions::eLinearAdvection);
647 DNekScalMat &AdvMat = *GetLocMatrix(advkey);
648
649 int rows = MassMat.GetRows();
650 int cols = MassMat.GetColumns();
651
652 DNekMatSharedPtr adr =
653 MemoryManager<DNekMat>::AllocateSharedPtr(rows, cols);
654
655 NekDouble one = 1.0;
656 (*adr) = -lambda * MassMat + AdvMat;
657
658 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, adr);
659
660 // Clear memory for time-dependent matrices
661 DropLocMatrix(advkey);
662 if (!massVarcoeffs.empty())
663 {
664 DropLocMatrix(masskey);
665 }
666 }
667 break;
668 case StdRegions::eLinearAdvectionDiffusionReaction:
669 {
670 NekDouble lambda = mkey.GetConstFactor(StdRegions::eFactorLambda);
671
672 // Construct mass matrix
673 // Check for mass-specific varcoeffs to avoid unncessary
674 // re-computation of the elemental matrix every time step
675 StdRegions::VarCoeffMap massVarcoeffs = StdRegions::NullVarCoeffMap;
676 if (mkey.HasVarCoeff(StdRegions::eVarCoeffMass))
677 {
678 massVarcoeffs[StdRegions::eVarCoeffMass] =
679 mkey.GetVarCoeff(StdRegions::eVarCoeffMass);
680 }
681 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this,
682 mkey.GetConstFactors(), massVarcoeffs);
683 DNekScalMat &MassMat = *GetLocMatrix(masskey);
684
685 // Construct laplacian matrix (Check for varcoeffs)
686 // Take all varcoeffs if one or more are detected
687 // TODO We might want to have a map
688 // from MatrixType to Vector of Varcoeffs and vice-versa
689 StdRegions::VarCoeffMap lapVarcoeffs = StdRegions::NullVarCoeffMap;
690 if ((mkey.HasVarCoeff(StdRegions::eVarCoeffLaplacian)) ||
691 (mkey.HasVarCoeff(StdRegions::eVarCoeffD00)) ||
692 (mkey.HasVarCoeff(StdRegions::eVarCoeffD01)) ||
693 (mkey.HasVarCoeff(StdRegions::eVarCoeffD10)) ||
694 (mkey.HasVarCoeff(StdRegions::eVarCoeffD02)) ||
695 (mkey.HasVarCoeff(StdRegions::eVarCoeffD20)) ||
696 (mkey.HasVarCoeff(StdRegions::eVarCoeffD11)) ||
697 (mkey.HasVarCoeff(StdRegions::eVarCoeffD12)) ||
698 (mkey.HasVarCoeff(StdRegions::eVarCoeffD21)) ||
699 (mkey.HasVarCoeff(StdRegions::eVarCoeffD22)))
700 {
701 lapVarcoeffs = mkey.GetVarCoeffs();
702 }
703 MatrixKey lapkey(StdRegions::eLaplacian, mkey.GetShapeType(), *this,
704 mkey.GetConstFactors(), lapVarcoeffs);
705 DNekScalMat &LapMat = *GetLocMatrix(lapkey);
706
707 // Construct advection matrix
708 // Check for varcoeffs not required;
709 // assume advection velocity is always time-dependent
710 MatrixKey advkey(mkey, StdRegions::eLinearAdvection);
711 DNekScalMat &AdvMat = *GetLocMatrix(advkey);
712
713 int rows = LapMat.GetRows();
714 int cols = LapMat.GetColumns();
715
716 DNekMatSharedPtr adr =
717 MemoryManager<DNekMat>::AllocateSharedPtr(rows, cols);
718
719 NekDouble one = 1.0;
720 (*adr) = LapMat - lambda * MassMat + AdvMat;
721
722 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, adr);
723
724 // Clear memory for time-dependent matrices
725 DropLocMatrix(advkey);
726 if (!massVarcoeffs.empty())
727 {
728 DropLocMatrix(masskey);
729 }
730 if (!lapVarcoeffs.empty())
731 {
732 DropLocMatrix(lapkey);
733 }
734 }
735 break;
736 case StdRegions::eLinearAdvectionDiffusionReactionGJP:
737 {
738 // Copied mostly from ADR solve to have fine-grain control
739 // over updating only advection matrix, relevant for performance!
740 NekDouble lambda = mkey.GetConstFactor(StdRegions::eFactorLambda);
741
742 // Construct mass matrix (Check for varcoeffs)
743 StdRegions::VarCoeffMap massVarcoeffs = StdRegions::NullVarCoeffMap;
744 if (mkey.HasVarCoeff(StdRegions::eVarCoeffMass))
745 {
746 massVarcoeffs[StdRegions::eVarCoeffMass] =
747 mkey.GetVarCoeff(StdRegions::eVarCoeffMass);
748 }
749 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this,
750 mkey.GetConstFactors(), massVarcoeffs);
751 DNekScalMat &MassMat = *GetLocMatrix(masskey);
752
753 // Construct laplacian matrix (Check for varcoeffs)
754 StdRegions::VarCoeffMap lapVarcoeffs = StdRegions::NullVarCoeffMap;
755 if ((mkey.HasVarCoeff(StdRegions::eVarCoeffLaplacian)) ||
756 (mkey.HasVarCoeff(StdRegions::eVarCoeffD00)) ||
757 (mkey.HasVarCoeff(StdRegions::eVarCoeffD01)) ||
758 (mkey.HasVarCoeff(StdRegions::eVarCoeffD10)) ||
759 (mkey.HasVarCoeff(StdRegions::eVarCoeffD02)) ||
760 (mkey.HasVarCoeff(StdRegions::eVarCoeffD20)) ||
761 (mkey.HasVarCoeff(StdRegions::eVarCoeffD11)) ||
762 (mkey.HasVarCoeff(StdRegions::eVarCoeffD12)) ||
763 (mkey.HasVarCoeff(StdRegions::eVarCoeffD21)) ||
764 (mkey.HasVarCoeff(StdRegions::eVarCoeffD22)))
765 {
766 lapVarcoeffs = mkey.GetVarCoeffs();
767 }
768 MatrixKey lapkey(StdRegions::eLaplacian, mkey.GetShapeType(), *this,
769 mkey.GetConstFactors(), lapVarcoeffs);
770 DNekScalMat &LapMat = *GetLocMatrix(lapkey);
771
772 // Construct advection matrix
773 // (assume advection velocity defined and non-zero)
774 // Could check L2(AdvectionVelocity) or HasVarCoeff
775 MatrixKey advkey(mkey, StdRegions::eLinearAdvection);
776 DNekScalMat &AdvMat = *GetLocMatrix(advkey);
777
778 // Generate a local copy of traceMat
779 MatrixKey gjpkey(StdRegions::eNormDerivOnTrace, mkey.GetShapeType(),
780 *this, mkey.GetConstFactors());
781 DNekScalMat &NDTraceMat = *GetLocMatrix(gjpkey);
782
783 NekDouble gjpfactor = mkey.GetConstFactor(StdRegions::eFactorGJP);
784 ASSERTL1(mkey.ConstFactorExists(StdRegions::eFactorGJP),
785 "Need to specify eFactorGJP to construct "
786 "a LinearAdvectionDiffusionReactionGJP matrix");
787
788 int rows = LapMat.GetRows();
789 int cols = LapMat.GetColumns();
790
791 DNekMatSharedPtr adr =
792 MemoryManager<DNekMat>::AllocateSharedPtr(rows, cols);
793
794 NekDouble one = 1.0;
795 (*adr) =
796 LapMat - lambda * MassMat + AdvMat + gjpfactor * NDTraceMat;
797
798 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, adr);
799
800 // Clear memory
801 DropLocMatrix(advkey);
802 DropLocMatrix(masskey);
803 DropLocMatrix(lapkey);
804 DropLocMatrix(gjpkey);
805 }
806 break;
807 case StdRegions::eNormDerivOnTrace:
808 {
809 NekDouble one = 1.0;
810 DNekMatSharedPtr mat = Expansion3D::v_GenMatrix(mkey);
811
812 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
813 }
814 break;
815 case StdRegions::eIProductWRTBase:
816 {
817 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
818 {
819 NekDouble one = 1.0;
820 DNekMatSharedPtr mat = GenMatrix(mkey);
821 returnval =
822 MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
823 }
824 else
825 {
826 NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
827 DNekMatSharedPtr mat = GetStdMatrix(mkey);
828 returnval =
829 MemoryManager<DNekScalMat>::AllocateSharedPtr(jac, mat);
830 }
831 }
832 break;
833 case StdRegions::eInvHybridDGHelmholtz:
834 {
835 NekDouble one = 1.0;
836
837 MatrixKey hkey(StdRegions::eHybridDGHelmholtz, DetShapeType(),
838 *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
839 DNekMatSharedPtr mat = GenMatrix(hkey);
840
841 mat->Invert();
842 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
843 }
844 break;
845 case StdRegions::ePreconLinearSpace:
846 {
847 NekDouble one = 1.0;
848 MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(),
849 *this, mkey.GetConstFactors(),
850 mkey.GetVarCoeffs());
851 DNekScalBlkMatSharedPtr helmStatCond =
852 GetLocStaticCondMatrix(helmkey);
853 DNekScalMatSharedPtr A = helmStatCond->GetBlock(0, 0);
854 DNekMatSharedPtr R = BuildVertexMatrix(A);
855
856 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, R);
857 }
858 break;
859 case StdRegions::ePreconLinearSpaceMass:
860 {
861 NekDouble one = 1.0;
862 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
863 DNekScalBlkMatSharedPtr massStatCond =
864 GetLocStaticCondMatrix(masskey);
865 DNekScalMatSharedPtr A = massStatCond->GetBlock(0, 0);
866 DNekMatSharedPtr R = BuildVertexMatrix(A);
867
868 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, R);
869 }
870 break;
871 case StdRegions::ePreconR:
872 {
873 NekDouble one = 1.0;
874 MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(),
875 *this, mkey.GetConstFactors(),
876 mkey.GetVarCoeffs());
877 DNekScalBlkMatSharedPtr helmStatCond =
878 GetLocStaticCondMatrix(helmkey);
879
880 DNekScalMatSharedPtr A = helmStatCond->GetBlock(0, 0);
881
882 DNekScalMatSharedPtr Atmp;
883
884 DNekMatSharedPtr R =
885 BuildTransformationMatrix(A, mkey.GetMatrixType());
886
887 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, R);
888
889 // Free memory,
890 // need to build keys to delete Mass and Laplacian
891 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
892 MatrixKey lapkey(StdRegions::eLaplacian, mkey.GetShapeType(), *this,
893 mkey.GetConstFactors(), mkey.GetVarCoeffs());
894 DropLocStaticCondMatrix(helmkey);
895 DropLocMatrix(helmkey);
896 DropLocMatrix(lapkey);
897 DropLocMatrix(masskey);
898 }
899 break;
900 case StdRegions::ePreconRMass:
901 {
902 NekDouble one = 1.0;
903 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
904 DNekScalBlkMatSharedPtr StatCond = GetLocStaticCondMatrix(masskey);
905 DNekScalMatSharedPtr A = StatCond->GetBlock(0, 0);
906
907 DNekScalMatSharedPtr Atmp;
908 DNekMatSharedPtr R =
909 BuildTransformationMatrix(A, mkey.GetMatrixType());
910
911 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, R);
912 }
913 break;
914 default:
915 {
916 NekDouble one = 1.0;
917 DNekMatSharedPtr mat = GenMatrix(mkey);
918
919 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
920 }
921 break;
922 }
923
924 return returnval;
925}
926
927/**
928 * Computes matrices needed for the HDG formulation. References to
929 * equations relate to the following paper (with a suitable changes in
930 * formulation to adapt to 3D):
931 * R. M. Kirby, S. J. Sherwin, B. Cockburn, To CG or to HDG: A
932 * Comparative Study, J. Sci. Comp P1-30
933 * DOI 10.1007/s10915-011-9501-7
934 * NOTE: VARIABLE COEFFICIENTS CASE IS NOT IMPLEMENTED
935 */
936DNekMatSharedPtr Expansion3D::v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
937{
938 // Variable coefficients are not implemented/////////
939 ASSERTL1(!mkey.HasVarCoeff(StdRegions::eVarCoeffD00),
940 "Matrix construction is not implemented for variable "
941 "coefficients at the moment");
942 ////////////////////////////////////////////////////
943
944 DNekMatSharedPtr returnval;
945
946 switch (mkey.GetMatrixType())
947 {
948 // (Z^e)^{-1} (Eqn. 33, P22)
949 case StdRegions::eHybridDGHelmholtz:
950 {
951 ASSERTL1(IsBoundaryInteriorExpansion(),
952 "HybridDGHelmholtz matrix not set up "
953 "for non boundary-interior expansions");
954
955 int i, j, k;
956 NekDouble lambdaval =
957 mkey.GetConstFactor(StdRegions::eFactorLambda);
958 NekDouble tau = mkey.GetConstFactor(StdRegions::eFactorTau);
959 int ncoeffs = GetNcoeffs();
960 int nfaces = GetNtraces();
961
962 Array<OneD, unsigned int> fmap;
963 Array<OneD, int> sign;
964 ExpansionSharedPtr FaceExp;
965 ExpansionSharedPtr FaceExp2;
966
967 int order_f, coordim = GetCoordim();
968 DNekScalMat &invMass = *GetLocMatrix(StdRegions::eInvMass);
969 StdRegions::MatrixType DerivType[3] = {StdRegions::eWeakDeriv0,
970 StdRegions::eWeakDeriv1,
971 StdRegions::eWeakDeriv2};
972
973 returnval =
974 MemoryManager<DNekMat>::AllocateSharedPtr(ncoeffs, ncoeffs);
975 DNekMat &Mat = *returnval;
976 Vmath::Zero(ncoeffs * ncoeffs, Mat.GetPtr(), 1);
977
978 // StdRegions::VarCoeffType Coeffs[3] = {StdRegions::eVarCoeffD00,
979 // StdRegions::eVarCoeffD11,
980 // StdRegions::eVarCoeffD22};
981 StdRegions::VarCoeffMap::const_iterator x;
982 const StdRegions::VarCoeffMap &varcoeffs = mkey.GetVarCoeffs();
983
984 for (i = 0; i < coordim; ++i)
985 {
986 if ((x = varcoeffs.find(StdRegions::eVarCoeffMF1x)) !=
987 varcoeffs.end())
988 {
989 StdRegions::VarCoeffMap VarCoeffDirDeriv;
990 VarCoeffDirDeriv[StdRegions::eVarCoeffMF] =
991 GetMF(i, coordim, varcoeffs);
992 VarCoeffDirDeriv[StdRegions::eVarCoeffMFDiv] =
993 GetMFDiv(i, varcoeffs);
994
995 MatrixKey Dmatkey(StdRegions::eWeakDirectionalDeriv,
996 DetShapeType(), *this,
997 StdRegions::NullConstFactorMap,
998 VarCoeffDirDeriv);
999
1000 DNekScalMat &Dmat = *GetLocMatrix(Dmatkey);
1001
1002 StdRegions::VarCoeffMap Weight;
1003 Weight[StdRegions::eVarCoeffMass] =
1004 GetMFMag(i, mkey.GetVarCoeffs());
1005
1006 MatrixKey invMasskey(StdRegions::eInvMass, DetShapeType(),
1007 *this, StdRegions::NullConstFactorMap,
1008 Weight);
1009
1010 DNekScalMat &invMass = *GetLocMatrix(invMasskey);
1011
1012 Mat = Mat + Dmat * invMass * Transpose(Dmat);
1013 }
1014 else
1015 {
1016 DNekScalMat &Dmat = *GetLocMatrix(DerivType[i]);
1017 Mat = Mat + Dmat * invMass * Transpose(Dmat);
1018 }
1019
1020 /*
1021 if(mkey.HasVarCoeff(Coeffs[i]))
1022 {
1023 MatrixKey DmatkeyL(DerivType[i], DetExpansionType(), *this,
1024 StdRegions::NullConstFactorMap,
1025 mkey.GetVarCoeffAsMap(Coeffs[i]));
1026 MatrixKey DmatkeyR(DerivType[i], DetExpansionType(), *this);
1027
1028 DNekScalMat &DmatL = *GetLocMatrix(DmatkeyL);
1029 DNekScalMat &DmatR = *GetLocMatrix(DmatkeyR);
1030 Mat = Mat + DmatL*invMass*Transpose(DmatR);
1031 }
1032 else
1033 {
1034 DNekScalMat &Dmat = *GetLocMatrix(DerivType[i]);
1035 Mat = Mat + Dmat*invMass*Transpose(Dmat);
1036 }
1037 */
1038 }
1039
1040 // Add Mass Matrix Contribution for Helmholtz problem
1041 DNekScalMat &Mass = *GetLocMatrix(StdRegions::eMass);
1042 Mat = Mat + lambdaval * Mass;
1043
1044 // Add tau*E_l using elemental mass matrices on each edge
1045 for (i = 0; i < nfaces; ++i)
1046 {
1047 FaceExp = GetTraceExp(i);
1048 order_f = FaceExp->GetNcoeffs();
1049
1050 IndexMapKey ikey(eFaceToElement, DetShapeType(),
1051 GetBasisNumModes(0), GetBasisNumModes(1),
1052 GetBasisNumModes(2), i, GetTraceOrient(i));
1053 IndexMapValuesSharedPtr map = GetIndexMap(ikey);
1054 ASSERTL0(order_f == (*map).size(),
1055 "HDGHelmholtz needs setting up for variable P");
1056
1057 // @TODO: Document
1058 /*
1059 StdRegions::VarCoeffMap edgeVarCoeffs;
1060 if (mkey.HasVarCoeff(StdRegions::eVarCoeffD00))
1061 {
1062 Array<OneD, NekDouble> mu(nq);
1063 GetPhysEdgeVarCoeffsFromElement(
1064 i, EdgeExp2,
1065 mkey.GetVarCoeff(StdRegions::eVarCoeffD00), mu);
1066 edgeVarCoeffs[StdRegions::eVarCoeffMass] = mu;
1067 }
1068 DNekScalMat &eMass = *EdgeExp->GetLocMatrix(
1069 StdRegions::eMass,
1070 StdRegions::NullConstFactorMap, edgeVarCoeffs);
1071 */
1072
1073 DNekScalMat &eMass = *FaceExp->GetLocMatrix(StdRegions::eMass);
1074
1075 for (j = 0; j < order_f; ++j)
1076 {
1077 for (k = 0; k < order_f; ++k)
1078 {
1079 Mat((*map)[j].index, (*map)[k].index) +=
1080 tau * (*map)[j].sign * (*map)[k].sign * eMass(j, k);
1081 }
1082 }
1083 }
1084 break;
1085 }
1086
1087 // U^e (P22)
1088 case StdRegions::eHybridDGLamToU:
1089 {
1090 int i, j, k;
1091 int nbndry = NumDGBndryCoeffs();
1092 int ncoeffs = GetNcoeffs();
1093 int nfaces = GetNtraces();
1094 NekDouble tau = mkey.GetConstFactor(StdRegions::eFactorTau);
1095
1096 Array<OneD, NekDouble> lambda(nbndry);
1097 DNekVec Lambda(nbndry, lambda, eWrapper);
1098 Array<OneD, NekDouble> ulam(ncoeffs);
1099 DNekVec Ulam(ncoeffs, ulam, eWrapper);
1100 Array<OneD, NekDouble> f(ncoeffs);
1101 DNekVec F(ncoeffs, f, eWrapper);
1102
1103 ExpansionSharedPtr FaceExp;
1104 // declare matrix space
1105 returnval =
1106 MemoryManager<DNekMat>::AllocateSharedPtr(ncoeffs, nbndry);
1107 DNekMat &Umat = *returnval;
1108
1109 // Z^e matrix
1110 MatrixKey newkey(StdRegions::eInvHybridDGHelmholtz, DetShapeType(),
1111 *this, mkey.GetConstFactors(),
1112 mkey.GetVarCoeffs());
1113 DNekScalMat &invHmat = *GetLocMatrix(newkey);
1114
1115 Array<OneD, unsigned int> fmap;
1116 Array<OneD, int> sign;
1117
1118 // alternative way to add boundary terms contribution
1119 int bndry_cnt = 0;
1120 for (i = 0; i < nfaces; ++i)
1121 {
1122 FaceExp = GetTraceExp(
1123 i); // temporary, need to rewrite AddHDGHelmholtzFaceTerms
1124 int nface = GetTraceNcoeffs(i);
1125 Array<OneD, NekDouble> face_lambda(nface);
1126
1127 const Array<OneD, const Array<OneD, NekDouble>> normals =
1128 GetTraceNormal(i);
1129
1130 for (j = 0; j < nface; ++j)
1131 {
1132 Vmath::Zero(nface, &face_lambda[0], 1);
1133 Vmath::Zero(ncoeffs, &f[0], 1);
1134 face_lambda[j] = 1.0;
1135
1136 SetFaceToGeomOrientation(i, face_lambda);
1137
1138 Array<OneD, NekDouble> tmp(FaceExp->GetTotPoints());
1139 FaceExp->BwdTrans(face_lambda, tmp);
1140 AddHDGHelmholtzFaceTerms(tau, i, tmp, mkey.GetVarCoeffs(),
1141 f);
1142
1143 Ulam = invHmat * F; // generate Ulam from lambda
1144
1145 // fill column of matrix
1146 for (k = 0; k < ncoeffs; ++k)
1147 {
1148 Umat(k, bndry_cnt) = Ulam[k];
1149 }
1150
1151 ++bndry_cnt;
1152 }
1153 }
1154
1155 //// Set up face expansions from local geom info
1156 // for(i = 0; i < nfaces; ++i)
1157 //{
1158 // FaceExp[i] = GetTraceExp(i);
1159 //}
1160 //
1161 //// for each degree of freedom of the lambda space
1162 //// calculate Umat entry
1163 //// Generate Lambda to U_lambda matrix
1164 // for(j = 0; j < nbndry; ++j)
1165 //{
1166 // // standard basis vectors e_j
1167 // Vmath::Zero(nbndry,&lambda[0],1);
1168 // Vmath::Zero(ncoeffs,&f[0],1);
1169 // lambda[j] = 1.0;
1170 //
1171 // //cout << Lambda;
1172 // SetTraceToGeomOrientation(lambda);
1173 // //cout << Lambda << endl;
1174 //
1175 // // Compute F = [I D_1 M^{-1} D_2 M^{-1}] C e_j
1176 // AddHDGHelmholtzTraceTerms(tau, lambda, FaceExp,
1177 // mkey.GetVarCoeffs(), f);
1178 //
1179 // // Compute U^e_j
1180 // Ulam = invHmat*F; // generate Ulam from lambda
1181 //
1182 // // fill column of matrix
1183 // for(k = 0; k < ncoeffs; ++k)
1184 // {
1185 // Umat(k,j) = Ulam[k];
1186 // }
1187 //}
1188 }
1189 break;
1190 // Q_0, Q_1, Q_2 matrices (P23)
1191 // Each are a product of a row of Eqn 32 with the C matrix.
1192 // Rather than explicitly computing all of Eqn 32, we note each
1193 // row is almost a multiple of U^e, so use that as our starting
1194 // point.
1195 case StdRegions::eHybridDGLamToQ0:
1196 case StdRegions::eHybridDGLamToQ1:
1197 case StdRegions::eHybridDGLamToQ2:
1198 {
1199 int i = 0;
1200 int j = 0;
1201 int k = 0;
1202 int dir = 0;
1203 int nbndry = NumDGBndryCoeffs();
1204 int coordim = GetCoordim();
1205 int ncoeffs = GetNcoeffs();
1206 int nfaces = GetNtraces();
1207
1208 Array<OneD, NekDouble> lambda(nbndry);
1209 DNekVec Lambda(nbndry, lambda, eWrapper);
1210 Array<OneD, ExpansionSharedPtr> FaceExp(nfaces);
1211
1212 Array<OneD, NekDouble> ulam(ncoeffs);
1213 DNekVec Ulam(ncoeffs, ulam, eWrapper);
1214 Array<OneD, NekDouble> f(ncoeffs);
1215 DNekVec F(ncoeffs, f, eWrapper);
1216
1217 // declare matrix space
1218 returnval =
1219 MemoryManager<DNekMat>::AllocateSharedPtr(ncoeffs, nbndry);
1220 DNekMat &Qmat = *returnval;
1221
1222 // Lambda to U matrix
1223 MatrixKey lamToUkey(StdRegions::eHybridDGLamToU, DetShapeType(),
1224 *this, mkey.GetConstFactors(),
1225 mkey.GetVarCoeffs());
1226 DNekScalMat &lamToU = *GetLocMatrix(lamToUkey);
1227
1228 // Inverse mass matrix
1229 DNekScalMat &invMass = *GetLocMatrix(StdRegions::eInvMass);
1230
1231 for (i = 0; i < nfaces; ++i)
1232 {
1233 FaceExp[i] = GetTraceExp(i);
1234 }
1235
1236 // Weak Derivative matrix
1237 DNekScalMatSharedPtr Dmat;
1238 switch (mkey.GetMatrixType())
1239 {
1240 case StdRegions::eHybridDGLamToQ0:
1241 dir = 0;
1242 Dmat = GetLocMatrix(StdRegions::eWeakDeriv0);
1243 break;
1244 case StdRegions::eHybridDGLamToQ1:
1245 dir = 1;
1246 Dmat = GetLocMatrix(StdRegions::eWeakDeriv1);
1247 break;
1248 case StdRegions::eHybridDGLamToQ2:
1249 dir = 2;
1250 Dmat = GetLocMatrix(StdRegions::eWeakDeriv2);
1251 break;
1252 default:
1253 ASSERTL0(false, "Direction not known");
1254 break;
1255 }
1256
1257 // DNekScalMatSharedPtr Dmat;
1258 // DNekScalMatSharedPtr &invMass;
1259 StdRegions::VarCoeffMap::const_iterator x;
1260 const StdRegions::VarCoeffMap &varcoeffs = mkey.GetVarCoeffs();
1261 if ((x = varcoeffs.find(StdRegions::eVarCoeffMF1x)) !=
1262 varcoeffs.end())
1263 {
1264 StdRegions::VarCoeffMap VarCoeffDirDeriv;
1265 VarCoeffDirDeriv[StdRegions::eVarCoeffMF] =
1266 GetMF(dir, coordim, varcoeffs);
1267 VarCoeffDirDeriv[StdRegions::eVarCoeffMFDiv] =
1268 GetMFDiv(dir, varcoeffs);
1269
1270 MatrixKey Dmatkey(
1271 StdRegions::eWeakDirectionalDeriv, DetShapeType(), *this,
1272 StdRegions::NullConstFactorMap, VarCoeffDirDeriv);
1273
1274 Dmat = GetLocMatrix(Dmatkey);
1275
1276 StdRegions::VarCoeffMap Weight;
1277 Weight[StdRegions::eVarCoeffMass] =
1278 GetMFMag(dir, mkey.GetVarCoeffs());
1279
1280 MatrixKey invMasskey(StdRegions::eInvMass, DetShapeType(),
1281 *this, StdRegions::NullConstFactorMap,
1282 Weight);
1283
1284 invMass = *GetLocMatrix(invMasskey);
1285 }
1286 else
1287 {
1288 // Inverse mass matrix
1289 invMass = *GetLocMatrix(StdRegions::eInvMass);
1290 }
1291
1292 // for each degree of freedom of the lambda space
1293 // calculate Qmat entry
1294 // Generate Lambda to Q_lambda matrix
1295 for (j = 0; j < nbndry; ++j)
1296 {
1297 Vmath::Zero(nbndry, &lambda[0], 1);
1298 lambda[j] = 1.0;
1299
1300 // for lambda[j] = 1 this is the solution to ulam
1301 for (k = 0; k < ncoeffs; ++k)
1302 {
1303 Ulam[k] = lamToU(k, j);
1304 }
1305
1306 // -D^T ulam
1307 Vmath::Neg(ncoeffs, &ulam[0], 1);
1308 F = Transpose(*Dmat) * Ulam;
1309
1310 SetTraceToGeomOrientation(lambda);
1311
1312 // Add the C terms resulting from the I's on the
1313 // diagonals of Eqn 32
1314 AddNormTraceInt(dir, lambda, FaceExp, f, mkey.GetVarCoeffs());
1315
1316 // finally multiply by inverse mass matrix
1317 Ulam = invMass * F;
1318
1319 // fill column of matrix (Qmat is in column major format)
1320 Vmath::Vcopy(ncoeffs, &ulam[0], 1,
1321 &(Qmat.GetPtr())[0] + j * ncoeffs, 1);
1322 }
1323 }
1324 break;
1325 // Matrix K (P23)
1326 case StdRegions::eHybridDGHelmBndLam:
1327 {
1328 int i, j, f, cnt;
1329 int order_f, nquad_f;
1330 int nbndry = NumDGBndryCoeffs();
1331 int nfaces = GetNtraces();
1332 NekDouble tau = mkey.GetConstFactor(StdRegions::eFactorTau);
1333
1334 Array<OneD, NekDouble> work, varcoeff_work;
1335 Array<OneD, const Array<OneD, NekDouble>> normals;
1336 Array<OneD, ExpansionSharedPtr> FaceExp(nfaces);
1337 Array<OneD, NekDouble> lam(nbndry);
1338
1339 Array<OneD, unsigned int> fmap;
1340 Array<OneD, int> sign;
1341
1342 // declare matrix space
1343 returnval =
1344 MemoryManager<DNekMat>::AllocateSharedPtr(nbndry, nbndry);
1345 DNekMat &BndMat = *returnval;
1346
1347 DNekScalMatSharedPtr LamToQ[3];
1348
1349 // Matrix to map Lambda to U
1350 MatrixKey LamToUkey(StdRegions::eHybridDGLamToU, DetShapeType(),
1351 *this, mkey.GetConstFactors(),
1352 mkey.GetVarCoeffs());
1353 DNekScalMat &LamToU = *GetLocMatrix(LamToUkey);
1354
1355 // Matrix to map Lambda to Q0
1356 MatrixKey LamToQ0key(StdRegions::eHybridDGLamToQ0, DetShapeType(),
1357 *this, mkey.GetConstFactors(),
1358 mkey.GetVarCoeffs());
1359 LamToQ[0] = GetLocMatrix(LamToQ0key);
1360
1361 // Matrix to map Lambda to Q1
1362 MatrixKey LamToQ1key(StdRegions::eHybridDGLamToQ1, DetShapeType(),
1363 *this, mkey.GetConstFactors(),
1364 mkey.GetVarCoeffs());
1365 LamToQ[1] = GetLocMatrix(LamToQ1key);
1366
1367 // Matrix to map Lambda to Q2
1368 MatrixKey LamToQ2key(StdRegions::eHybridDGLamToQ2, DetShapeType(),
1369 *this, mkey.GetConstFactors(),
1370 mkey.GetVarCoeffs());
1371 LamToQ[2] = GetLocMatrix(LamToQ2key);
1372
1373 // Set up edge segment expansions from local geom info
1374 const StdRegions::VarCoeffMap &varcoeffs = mkey.GetVarCoeffs();
1375 for (i = 0; i < nfaces; ++i)
1376 {
1377 FaceExp[i] = GetTraceExp(i);
1378 }
1379
1380 // Set up matrix derived from <mu, Q_lam.n - \tau (U_lam - Lam) >
1381 for (i = 0; i < nbndry; ++i)
1382 {
1383 cnt = 0;
1384
1385 Vmath::Zero(nbndry, lam, 1);
1386 lam[i] = 1.0;
1387 SetTraceToGeomOrientation(lam);
1388
1389 for (f = 0; f < nfaces; ++f)
1390 {
1391 order_f = FaceExp[f]->GetNcoeffs();
1392 nquad_f = FaceExp[f]->GetNumPoints(0) *
1393 FaceExp[f]->GetNumPoints(1);
1394 normals = GetTraceNormal(f);
1395
1396 work = Array<OneD, NekDouble>(nquad_f);
1397 varcoeff_work = Array<OneD, NekDouble>(nquad_f);
1398
1399 IndexMapKey ikey(eFaceToElement, DetShapeType(),
1400 GetBasisNumModes(0), GetBasisNumModes(1),
1401 GetBasisNumModes(2), f, GetTraceOrient(f));
1402 IndexMapValuesSharedPtr map = GetIndexMap(ikey);
1403
1404 // @TODO Variable coefficients
1405 /*
1406 StdRegions::VarCoeffType VarCoeff[3] =
1407 {StdRegions::eVarCoeffD00, StdRegions::eVarCoeffD11,
1408 StdRegions::eVarCoeffD22};
1409 const StdRegions::VarCoeffMap &varcoeffs =
1410 mkey.GetVarCoeffs(); StdRegions::VarCoeffMap::const_iterator
1411 x;
1412 */
1413
1414 // Q0 * n0 (BQ_0 terms)
1415 Array<OneD, NekDouble> faceCoeffs(order_f);
1416 Array<OneD, NekDouble> facePhys(nquad_f);
1417 for (j = 0; j < order_f; ++j)
1418 {
1419 faceCoeffs[j] =
1420 (*map)[j].sign * (*LamToQ[0])((*map)[j].index, i);
1421 }
1422
1423 FaceExp[f]->BwdTrans(faceCoeffs, facePhys);
1424
1425 // @TODO Variable coefficients
1426 // Multiply by variable coefficient
1427 /*
1428 if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
1429 {
1430 GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
1431 Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
1432 }
1433 */
1434
1435 if (varcoeffs.find(StdRegions::eVarCoeffMF1x) !=
1436 varcoeffs.end())
1437 {
1438 Array<OneD, NekDouble> ncdotMF = GetnFacecdotMF(
1439 0, f, FaceExp[f], normals, varcoeffs);
1440
1441 Vmath::Vmul(nquad_f, ncdotMF, 1, facePhys, 1, work, 1);
1442 }
1443 else
1444 {
1445 Vmath::Vmul(nquad_f, normals[0], 1, facePhys, 1, work,
1446 1);
1447 }
1448
1449 // Q1 * n1 (BQ_1 terms)
1450 for (j = 0; j < order_f; ++j)
1451 {
1452 faceCoeffs[j] =
1453 (*map)[j].sign * (*LamToQ[1])((*map)[j].index, i);
1454 }
1455
1456 FaceExp[f]->BwdTrans(faceCoeffs, facePhys);
1457
1458 // @TODO Variable coefficients
1459 // Multiply by variable coefficients
1460 /*
1461 if ((x = varcoeffs.find(VarCoeff[1])) != varcoeffs.end())
1462 {
1463 GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
1464 Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
1465 }
1466 */
1467
1468 if ((varcoeffs.find(StdRegions::eVarCoeffMF1x)) !=
1469 varcoeffs.end())
1470 {
1471 Array<OneD, NekDouble> ncdotMF = GetnFacecdotMF(
1472 1, f, FaceExp[f], normals, varcoeffs);
1473
1474 Vmath::Vvtvp(nquad_f, ncdotMF, 1, facePhys, 1, work, 1,
1475 work, 1);
1476 }
1477 else
1478 {
1479 Vmath::Vvtvp(nquad_f, normals[1], 1, facePhys, 1, work,
1480 1, work, 1);
1481 }
1482
1483 // Q2 * n2 (BQ_2 terms)
1484 for (j = 0; j < order_f; ++j)
1485 {
1486 faceCoeffs[j] =
1487 (*map)[j].sign * (*LamToQ[2])((*map)[j].index, i);
1488 }
1489
1490 FaceExp[f]->BwdTrans(faceCoeffs, facePhys);
1491
1492 // @TODO Variable coefficients
1493 // Multiply by variable coefficients
1494 /*
1495 if ((x = varcoeffs.find(VarCoeff[2])) != varcoeffs.end())
1496 {
1497 GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
1498 Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
1499 }
1500 */
1501
1502 if (varcoeffs.find(StdRegions::eVarCoeffMF1x) !=
1503 varcoeffs.end())
1504 {
1505 Array<OneD, NekDouble> ncdotMF = GetnFacecdotMF(
1506 2, f, FaceExp[f], normals, varcoeffs);
1507
1508 Vmath::Vvtvp(nquad_f, ncdotMF, 1, facePhys, 1, work, 1,
1509 work, 1);
1510 }
1511 else
1512 {
1513 Vmath::Vvtvp(nquad_f, normals[2], 1, facePhys, 1, work,
1514 1, work, 1);
1515 }
1516
1517 // - tau (ulam - lam)
1518 // Corresponds to the G and BU terms.
1519 for (j = 0; j < order_f; ++j)
1520 {
1521 faceCoeffs[j] =
1522 (*map)[j].sign * LamToU((*map)[j].index, i) -
1523 lam[cnt + j];
1524 }
1525
1526 FaceExp[f]->BwdTrans(faceCoeffs, facePhys);
1527
1528 // @TODO Variable coefficients
1529 // Multiply by variable coefficients
1530 /*
1531 if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
1532 {
1533 GetPhysEdgeVarCoeffsFromElement(e,FaceExp[f],x->second,varcoeff_work);
1534 Vmath::Vmul(nquad_f,varcoeff_work,1,FaceExp[f]->GetPhys(),1,FaceExp[f]->UpdatePhys(),1);
1535 }
1536 */
1537
1538 Vmath::Svtvp(nquad_f, -tau, facePhys, 1, work, 1, work, 1);
1539
1540 // @TODO Add variable coefficients
1541 FaceExp[f]->IProductWRTBase(work, faceCoeffs);
1542
1543 SetFaceToGeomOrientation(f, faceCoeffs);
1544
1545 for (j = 0; j < order_f; ++j)
1546 {
1547 BndMat(cnt + j, i) = faceCoeffs[j];
1548 }
1549
1550 cnt += order_f;
1551 }
1552 }
1553 }
1554 break;
1555 // HDG postprocessing
1556 case StdRegions::eInvLaplacianWithUnityMean:
1557 {
1558 MatrixKey lapkey(StdRegions::eLaplacian, DetShapeType(), *this,
1559 mkey.GetConstFactors(), mkey.GetVarCoeffs());
1560 DNekScalMat &LapMat = *GetLocMatrix(lapkey);
1561
1562 returnval = MemoryManager<DNekMat>::AllocateSharedPtr(
1563 LapMat.GetRows(), LapMat.GetColumns());
1564 DNekMatSharedPtr lmat = returnval;
1565
1566 (*lmat) = LapMat;
1567
1568 // replace first column with inner product wrt 1
1569 int nq = GetTotPoints();
1570 Array<OneD, NekDouble> tmp(nq);
1571 Array<OneD, NekDouble> outarray(m_ncoeffs);
1572 Vmath::Fill(nq, 1.0, tmp, 1);
1573 IProductWRTBase(tmp, outarray);
1574
1575 Vmath::Vcopy(m_ncoeffs, &outarray[0], 1, &(lmat->GetPtr())[0], 1);
1576
1577 lmat->Invert();
1578 }
1579 break;
1580 case StdRegions::eNormDerivOnTrace:
1581 {
1582 int ntraces = GetNtraces();
1583 int ncoords = GetCoordim();
1584 int nphys = GetTotPoints();
1585 Array<OneD, const Array<OneD, NekDouble>> normals;
1586 Array<OneD, NekDouble> phys(nphys);
1587 returnval =
1588 MemoryManager<DNekMat>::AllocateSharedPtr(m_ncoeffs, m_ncoeffs);
1589 DNekMat &Mat = *returnval;
1590 Vmath::Zero(m_ncoeffs * m_ncoeffs, Mat.GetPtr(), 1);
1591
1592 Array<OneD, Array<OneD, NekDouble>> Deriv(3, NullNekDouble1DArray);
1593
1594 for (int d = 0; d < ncoords; ++d)
1595 {
1596 Deriv[d] = Array<OneD, NekDouble>(nphys);
1597 }
1598
1599 Array<OneD, int> tracepts(ntraces);
1600 Array<OneD, ExpansionSharedPtr> traceExp(ntraces);
1601 int maxtpts = 0;
1602 for (int t = 0; t < ntraces; ++t)
1603 {
1604 traceExp[t] = GetTraceExp(t);
1605 tracepts[t] = traceExp[t]->GetTotPoints();
1606 maxtpts = (maxtpts > tracepts[t]) ? maxtpts : tracepts[t];
1607 }
1608
1609 Array<OneD, NekDouble> val(maxtpts), tmp, tmp1;
1610
1611 Array<OneD, Array<OneD, NekDouble>> dphidn(ntraces);
1612 for (int t = 0; t < ntraces; ++t)
1613 {
1614 dphidn[t] =
1615 Array<OneD, NekDouble>(m_ncoeffs * tracepts[t], 0.0);
1616 }
1617
1618 for (int i = 0; i < m_ncoeffs; ++i)
1619 {
1620 FillMode(i, phys);
1621 PhysDeriv(phys, Deriv[0], Deriv[1], Deriv[2]);
1622
1623 for (int t = 0; t < ntraces; ++t)
1624 {
1625 const NormalVector norm = GetTraceNormal(t);
1626
1627 LibUtilities::BasisKey fromkey0 = GetTraceBasisKey(t, 0);
1628 LibUtilities::BasisKey fromkey1 = GetTraceBasisKey(t, 1);
1629 LibUtilities::BasisKey tokey0 =
1630 traceExp[t]->GetBasis(0)->GetBasisKey();
1631 LibUtilities::BasisKey tokey1 =
1632 traceExp[t]->GetBasis(1)->GetBasisKey();
1633 bool DoInterp =
1634 (fromkey0 != tokey0) || (fromkey1 != tokey1);
1635
1636 // for variable p need add check and
1637 // interpolation here.
1638
1639 Array<OneD, NekDouble> n(tracepts[t]);
1640 ;
1641 for (int d = 0; d < ncoords; ++d)
1642 {
1643 // if variable p may need to interpolate
1644 if (DoInterp)
1645 {
1646 LibUtilities::Interp2D(fromkey0, fromkey1, norm[d],
1647 tokey0, tokey1, n);
1648 }
1649 else
1650 {
1651 n = norm[d];
1652 }
1653
1654 GetTracePhysVals(t, traceExp[t], Deriv[d], val,
1655 v_GetTraceOrient(t));
1656
1657 Vmath::Vvtvp(tracepts[t], n, 1, val, 1,
1658 tmp = dphidn[t] + i * tracepts[t], 1,
1659 tmp1 = dphidn[t] + i * tracepts[t], 1);
1660 }
1661 }
1662 }
1663
1664 for (int t = 0; t < ntraces; ++t)
1665 {
1666 int nt = tracepts[t];
1667 NekDouble h, p;
1668 TraceNormLen(t, h, p);
1669
1670 // scaling from GJP paper
1671 NekDouble scale =
1672 (p == 1) ? 0.02 * h * h : 0.8 * pow(p + 1, -4.0) * h * h;
1673
1674 for (int i = 0; i < m_ncoeffs; ++i)
1675 {
1676 for (int j = i; j < m_ncoeffs; ++j)
1677 {
1678 Vmath::Vmul(nt, dphidn[t] + i * nt, 1,
1679 dphidn[t] + j * nt, 1, val, 1);
1680 Mat(i, j) =
1681 Mat(i, j) + scale * traceExp[t]->Integral(val);
1682 }
1683 }
1684 }
1685
1686 // fill in symmetric components.
1687 for (int i = 0; i < m_ncoeffs; ++i)
1688 {
1689 for (int j = 0; j < i; ++j)
1690 {
1691 Mat(i, j) = Mat(j, i);
1692 }
1693 }
1694 }
1695 break;
1696 default:
1697 ASSERTL0(false,
1698 "This matrix type cannot be generated from this class");
1699 break;
1700 }
1701
1702 return returnval;
1703}
1704
1705void Expansion3D::v_AddFaceNormBoundaryInt(
1706 const int face, const ExpansionSharedPtr &FaceExp,
1707 const Array<OneD, const NekDouble> &Fn, Array<OneD, NekDouble> &outarray)
1708{
1709 int i, j;
1710
1711 /*
1712 * Coming into this routine, the velocity V will have been
1713 * multiplied by the trace normals to give the input vector Vn. By
1714 * convention, these normals are inwards facing for elements which
1715 * have FaceExp as their right-adjacent face. This conditional
1716 * statement therefore determines whether the normals must be
1717 * negated, since the integral being performed here requires an
1718 * outwards facing normal.
1719 */
1720 if (m_requireNeg.size() == 0)
1721 {
1722 m_requireNeg.resize(GetNtraces());
1723
1724 for (i = 0; i < GetNtraces(); ++i)
1725 {
1726 m_requireNeg[i] = false;
1727
1728 ExpansionSharedPtr faceExp = m_traceExp[i].lock();
1729
1730 if (faceExp->GetRightAdjacentElementExp())
1731 {
1732 if (faceExp->GetRightAdjacentElementExp()
1733 ->GetGeom()
1734 ->GetGlobalID() == GetGeom()->GetGlobalID())
1735 {
1736 m_requireNeg[i] = true;
1737 }
1738 }
1739 }
1740 }
1741
1742 IndexMapKey ikey(eFaceToElement, DetShapeType(), GetBasisNumModes(0),
1743 GetBasisNumModes(1), GetBasisNumModes(2), face,
1744 GetTraceOrient(face));
1745 IndexMapValuesSharedPtr map = GetIndexMap(ikey);
1746
1747 int order_e = (*map).size(); // Order of the element
1748 int n_coeffs = FaceExp->GetNcoeffs();
1749
1750 Array<OneD, NekDouble> faceCoeffs(n_coeffs);
1751
1752 if (n_coeffs != order_e) // Going to orthogonal space
1753 {
1754 FaceExp->FwdTrans(Fn, faceCoeffs);
1755
1756 int NumModesElementMax = FaceExp->GetBasis(0)->GetNumModes();
1757 int NumModesElementMin = m_base[0]->GetNumModes();
1758
1759 FaceExp->ReduceOrderCoeffs(NumModesElementMin, faceCoeffs, faceCoeffs);
1760
1761 StdRegions::StdMatrixKey masskey(StdRegions::eMass,
1762 FaceExp->DetShapeType(), *FaceExp);
1763 FaceExp->MassMatrixOp(faceCoeffs, faceCoeffs, masskey);
1764
1765 // Reorder coefficients for the lower degree face.
1766 int offset1 = 0, offset2 = 0;
1767
1768 if (FaceExp->DetShapeType() == LibUtilities::eQuadrilateral)
1769 {
1770 for (i = 0; i < NumModesElementMin; ++i)
1771 {
1772 for (j = 0; j < NumModesElementMin; ++j)
1773 {
1774 faceCoeffs[offset1 + j] = faceCoeffs[offset2 + j];
1775 }
1776 offset1 += NumModesElementMin;
1777 offset2 += NumModesElementMax;
1778 }
1779
1780 // Extract lower degree modes. TODO: Check this is correct.
1781 for (i = NumModesElementMin; i < NumModesElementMax; ++i)
1782 {
1783 for (j = NumModesElementMin; j < NumModesElementMax; ++j)
1784 {
1785 faceCoeffs[i * NumModesElementMax + j] = 0.0;
1786 }
1787 }
1788 }
1789
1790 if (FaceExp->DetShapeType() == LibUtilities::eTriangle)
1791 {
1792
1793 // Reorder coefficients for the lower degree face.
1794 int offset1 = 0, offset2 = 0;
1795
1796 for (i = 0; i < NumModesElementMin; ++i)
1797 {
1798 for (j = 0; j < NumModesElementMin - i; ++j)
1799 {
1800 faceCoeffs[offset1 + j] = faceCoeffs[offset2 + j];
1801 }
1802 offset1 += NumModesElementMin - i;
1803 offset2 += NumModesElementMax - i;
1804 }
1805 }
1806 }
1807 else
1808 {
1809 FaceExp->IProductWRTBase(Fn, faceCoeffs);
1810 }
1811
1812 if (m_requireNeg[face])
1813 {
1814 for (i = 0; i < order_e; ++i)
1815 {
1816 outarray[(*map)[i].index] -= (*map)[i].sign * faceCoeffs[i];
1817 }
1818 }
1819 else
1820 {
1821 for (i = 0; i < order_e; ++i)
1822 {
1823 outarray[(*map)[i].index] += (*map)[i].sign * faceCoeffs[i];
1824 }
1825 }
1826}
1827
1828/**
1829 * @brief Evaluate coefficients of weak deriviative in the direction dir
1830 * given the input coefficicents incoeffs and the imposed boundary
1831 * values in EdgeExp (which will have its phys space updated).
1832 */
1833void Expansion3D::v_DGDeriv(int dir,
1834 const Array<OneD, const NekDouble> &incoeffs,
1835 Array<OneD, ExpansionSharedPtr> &FaceExp,
1836 Array<OneD, Array<OneD, NekDouble>> &faceCoeffs,
1837 Array<OneD, NekDouble> &out_d)
1838{
1839 int ncoeffs = GetNcoeffs();
1840 StdRegions::MatrixType DerivType[3] = {StdRegions::eWeakDeriv0,
1841 StdRegions::eWeakDeriv1,
1842 StdRegions::eWeakDeriv2};
1843
1844 DNekScalMat &InvMass = *GetLocMatrix(StdRegions::eInvMass);
1845 DNekScalMat &Dmat = *GetLocMatrix(DerivType[dir]);
1846
1847 Array<OneD, NekDouble> coeffs = incoeffs;
1848 DNekVec Coeffs(ncoeffs, coeffs, eWrapper);
1849
1850 Coeffs = Transpose(Dmat) * Coeffs;
1851 Vmath::Neg(ncoeffs, coeffs, 1);
1852
1853 // Add the boundary integral including the relevant part of
1854 // the normal
1855 AddNormTraceInt(dir, FaceExp, faceCoeffs, coeffs);
1856
1857 DNekVec Out_d(ncoeffs, out_d, eWrapper);
1858
1859 Out_d = InvMass * Coeffs;
1860}
1861
1862void Expansion3D::v_AddRobinMassMatrix(
1863 const int face, const Array<OneD, const NekDouble> &primCoeffs,
1864 DNekMatSharedPtr &inoutmat)
1865{
1866 ASSERTL1(IsBoundaryInteriorExpansion(),
1867 "Not set up for non boundary-interior expansions");
1868 ASSERTL1(inoutmat->GetRows() == inoutmat->GetColumns(),
1869 "Assuming that input matrix was square");
1870
1871 int i, j;
1872 int id1, id2;
1873 ExpansionSharedPtr faceExp = m_traceExp[face].lock();
1874 int order_f = faceExp->GetNcoeffs();
1875
1876 Array<OneD, unsigned int> map;
1877 Array<OneD, int> sign;
1878
1879 StdRegions::VarCoeffMap varcoeffs;
1880 varcoeffs[StdRegions::eVarCoeffMass] = primCoeffs;
1881
1882 LibUtilities::ShapeType shapeType = faceExp->DetShapeType();
1883
1884 LocalRegions::MatrixKey mkey(StdRegions::eMass, shapeType, *faceExp,
1885 StdRegions::NullConstFactorMap, varcoeffs);
1886
1887 DNekScalMat &facemat = *faceExp->GetLocMatrix(mkey);
1888
1889 // Now need to identify a map which takes the local face
1890 // mass matrix to the matrix stored in inoutmat;
1891 // This can currently be deduced from the size of the matrix
1892
1893 // - if inoutmat.m_rows() == v_NCoeffs() it is a full
1894 // matrix system
1895
1896 // - if inoutmat.m_rows() == v_GetNverts() it is a vertex space
1897 // preconditioner.
1898
1899 // - if inoutmat.m_rows() == v_NumBndCoeffs() it is a
1900 // boundary CG system
1901
1902 // - if inoutmat.m_rows() == v_NumDGBndCoeffs() it is a
1903 // trace DG system; still needs implementing.
1904 int rows = inoutmat->GetRows();
1905
1906 if (rows == GetNcoeffs())
1907 {
1908 GetTraceToElementMap(face, map, sign, GetTraceOrient(face));
1909 }
1910 else if (rows == GetNverts())
1911 {
1912 int nfvert = faceExp->GetNverts();
1913
1914 // Need to find where linear vertices are in facemat
1915 Array<OneD, unsigned int> linmap;
1916 Array<OneD, int> linsign;
1917
1918 // Use a linear expansion to get correct mapping
1919 GetLinStdExp()->GetTraceToElementMap(face, linmap, linsign,
1920 GetTraceOrient(face));
1921
1922 // zero out sign map to remove all other modes
1923 sign = Array<OneD, int>(order_f, 0);
1924 map = Array<OneD, unsigned int>(order_f, (unsigned int)0);
1925
1926 int fmap;
1927 // Reset sign map to only have contribution from vertices
1928 for (i = 0; i < nfvert; ++i)
1929 {
1930 fmap = faceExp->GetVertexMap(i, true);
1931 sign[fmap] = 1;
1932
1933 // need to reset map
1934 map[fmap] = linmap[i];
1935 }
1936 }
1937 else if (rows == NumBndryCoeffs())
1938 {
1939 int nbndry = NumBndryCoeffs();
1940 Array<OneD, unsigned int> bmap(nbndry);
1941
1942 GetTraceToElementMap(face, map, sign, GetTraceOrient(face));
1943 GetBoundaryMap(bmap);
1944
1945 for (i = 0; i < order_f; ++i)
1946 {
1947 for (j = 0; j < nbndry; ++j)
1948 {
1949 if (map[i] == bmap[j])
1950 {
1951 map[i] = j;
1952 break;
1953 }
1954 }
1955 ASSERTL1(j != nbndry, "Did not find number in map");
1956 }
1957 }
1958 else if (rows == NumDGBndryCoeffs())
1959 {
1960 // possibly this should be a separate method
1961 int cnt = 0;
1962 map = Array<OneD, unsigned int>(order_f);
1963 sign = Array<OneD, int>(order_f, 1);
1964
1965 IndexMapKey ikey1(eFaceToElement, DetShapeType(), GetBasisNumModes(0),
1966 GetBasisNumModes(1), GetBasisNumModes(2), face,
1967 GetTraceOrient(face));
1968 IndexMapValuesSharedPtr map1 = GetIndexMap(ikey1);
1969 IndexMapKey ikey2(eFaceToElement, DetShapeType(), GetBasisNumModes(0),
1970 GetBasisNumModes(1), GetBasisNumModes(2), face,
1971 StdRegions::eDir1FwdDir1_Dir2FwdDir2);
1972 IndexMapValuesSharedPtr map2 = GetIndexMap(ikey2);
1973
1974 ASSERTL1((*map1).size() == (*map2).size(),
1975 "There is an error with the GetTraceToElementMap");
1976
1977 for (i = 0; i < face; ++i)
1978 {
1979 cnt += GetTraceNcoeffs(i);
1980 }
1981
1982 for (i = 0; i < (*map1).size(); ++i)
1983 {
1984 int idx = -1;
1985
1986 for (j = 0; j < (*map2).size(); ++j)
1987 {
1988 if ((*map1)[i].index == (*map2)[j].index)
1989 {
1990 idx = j;
1991 break;
1992 }
1993 }
1994
1995 ASSERTL2(idx >= 0, "Index not found");
1996 map[i] = idx + cnt;
1997 sign[i] = (*map2)[idx].sign;
1998 }
1999 }
2000 else
2001 {
2002 ASSERTL0(false, "Could not identify matrix type from dimension");
2003 }
2004
2005 for (i = 0; i < order_f; ++i)
2006 {
2007 id1 = map[i];
2008 for (j = 0; j < order_f; ++j)
2009 {
2010 id2 = map[j];
2011 (*inoutmat)(id1, id2) += facemat(i, j) * sign[i] * sign[j];
2012 }
2013 }
2014}
2015
2016DNekMatSharedPtr Expansion3D::v_BuildVertexMatrix(
2017 const DNekScalMatSharedPtr &r_bnd)
2018{
2019 MatrixStorage storage = eFULL;
2020 DNekMatSharedPtr vertexmatrix;
2021
2022 int nVerts, vid1, vid2, vMap1, vMap2;
2023 NekDouble VertexValue;
2024
2025 nVerts = GetNverts();
2026
2027 vertexmatrix =
2028 MemoryManager<DNekMat>::AllocateSharedPtr(nVerts, nVerts, 0.0, storage);
2029 DNekMat &VertexMat = (*vertexmatrix);
2030
2031 for (vid1 = 0; vid1 < nVerts; ++vid1)
2032 {
2033 vMap1 = GetVertexMap(vid1, true);
2034
2035 for (vid2 = 0; vid2 < nVerts; ++vid2)
2036 {
2037 vMap2 = GetVertexMap(vid2, true);
2038 VertexValue = (*r_bnd)(vMap1, vMap2);
2039 VertexMat.SetValue(vid1, vid2, VertexValue);
2040 }
2041 }
2042
2043 return vertexmatrix;
2044}
2045
2046DNekMatSharedPtr Expansion3D::v_BuildTransformationMatrix(
2047 const DNekScalMatSharedPtr &r_bnd, const StdRegions::MatrixType matrixType)
2048{
2049 int nVerts, nEdges;
2050 int eid, fid, vid, n, i;
2051
2052 int nBndCoeffs = NumBndryCoeffs();
2053
2054 const SpatialDomains::Geometry3D *geom = GetGeom3D();
2055
2056 // Get geometric information about this element
2057 nVerts = GetNverts();
2058 nEdges = GetNedges();
2059
2060 /*************************************/
2061 /* Vetex-edge & vertex-face matrices */
2062 /*************************************/
2063
2064 /**
2065 * The matrix component of \f$\mathbf{R}\f$ is given by \f[
2066 * \mathbf{R^{T}_{v}}=
2067 * -\mathbf{S}^{-1}_{ef,ef}\mathbf{S}^{T}_{v,ef}\f]
2068 *
2069 * For every vertex mode we extract the submatrices from statically
2070 * condensed matrix \f$\mathbf{S}\f$ corresponding to the coupling
2071 * between the attached edges and faces of a vertex
2072 * (\f$\mathbf{S_{ef,ef}}\f$). This matrix is then inverted and
2073 * multiplied by the submatrix representing the coupling between a
2074 * vertex and the attached edges and faces
2075 * (\f$\mathbf{S_{v,ef}}\f$).
2076 */
2077
2078 int nmodes;
2079 int m;
2080 NekDouble VertexEdgeFaceValue;
2081
2082 // The number of connected edges/faces is 3 (for all elements)
2083 int nConnectedEdges = 3;
2084 int nConnectedFaces = 3;
2085
2086 // Location in the matrix
2087 Array<OneD, Array<OneD, unsigned int>> MatEdgeLocation(nConnectedEdges);
2088 Array<OneD, Array<OneD, unsigned int>> MatFaceLocation(nConnectedFaces);
2089
2090 // Define storage for vertex transpose matrix and zero all entries
2091 MatrixStorage storage = eFULL;
2092 DNekMatSharedPtr transformationmatrix;
2093
2094 transformationmatrix = MemoryManager<DNekMat>::AllocateSharedPtr(
2095 nBndCoeffs, nBndCoeffs, 0.0, storage);
2096
2097 DNekMat &R = (*transformationmatrix);
2098
2099 // Build the vertex-edge/face transform matrix: This matrix is
2100 // constructed from the submatrices corresponding to the couping
2101 // between each vertex and the attached edges/faces
2102 for (vid = 0; vid < nVerts; ++vid)
2103 {
2104 // Row and column size of the vertex-edge/face matrix
2105 int efRow = GetEdgeNcoeffs(geom->GetVertexEdgeMap(vid, 0)) +
2106 GetEdgeNcoeffs(geom->GetVertexEdgeMap(vid, 1)) +
2107 GetEdgeNcoeffs(geom->GetVertexEdgeMap(vid, 2)) +
2108 GetTraceIntNcoeffs(geom->GetVertexFaceMap(vid, 0)) +
2109 GetTraceIntNcoeffs(geom->GetVertexFaceMap(vid, 1)) +
2110 GetTraceIntNcoeffs(geom->GetVertexFaceMap(vid, 2)) - 6;
2111
2112 int nedgemodesconnected =
2113 GetEdgeNcoeffs(geom->GetVertexEdgeMap(vid, 0)) +
2114 GetEdgeNcoeffs(geom->GetVertexEdgeMap(vid, 1)) +
2115 GetEdgeNcoeffs(geom->GetVertexEdgeMap(vid, 2)) - 6;
2116 Array<OneD, unsigned int> edgemodearray(nedgemodesconnected);
2117
2118 int nfacemodesconnected =
2119 GetTraceIntNcoeffs(geom->GetVertexFaceMap(vid, 0)) +
2120 GetTraceIntNcoeffs(geom->GetVertexFaceMap(vid, 1)) +
2121 GetTraceIntNcoeffs(geom->GetVertexFaceMap(vid, 2));
2122 Array<OneD, unsigned int> facemodearray(nfacemodesconnected);
2123
2124 int offset = 0;
2125
2126 // Create array of edge modes
2127 for (eid = 0; eid < nConnectedEdges; ++eid)
2128 {
2129 MatEdgeLocation[eid] =
2130 GetEdgeInverseBoundaryMap(geom->GetVertexEdgeMap(vid, eid));
2131 nmodes = MatEdgeLocation[eid].size();
2132
2133 if (nmodes)
2134 {
2135 Vmath::Vcopy(nmodes, &MatEdgeLocation[eid][0], 1,
2136 &edgemodearray[offset], 1);
2137 }
2138
2139 offset += nmodes;
2140 }
2141
2142 offset = 0;
2143
2144 // Create array of face modes
2145 for (fid = 0; fid < nConnectedFaces; ++fid)
2146 {
2147 MatFaceLocation[fid] =
2148 GetTraceInverseBoundaryMap(geom->GetVertexFaceMap(vid, fid));
2149 nmodes = MatFaceLocation[fid].size();
2150
2151 if (nmodes)
2152 {
2153 Vmath::Vcopy(nmodes, &MatFaceLocation[fid][0], 1,
2154 &facemodearray[offset], 1);
2155 }
2156 offset += nmodes;
2157 }
2158
2159 DNekMatSharedPtr vertexedgefacetransformmatrix =
2160 MemoryManager<DNekMat>::AllocateSharedPtr(1, efRow, 0.0, storage);
2161 DNekMat &Sveft = (*vertexedgefacetransformmatrix);
2162
2163 DNekMatSharedPtr vertexedgefacecoupling =
2164 MemoryManager<DNekMat>::AllocateSharedPtr(1, efRow, 0.0, storage);
2165 DNekMat &Svef = (*vertexedgefacecoupling);
2166
2167 // Vertex-edge coupling
2168 for (n = 0; n < nedgemodesconnected; ++n)
2169 {
2170 // Matrix value for each coefficient location
2171 VertexEdgeFaceValue = (*r_bnd)(GetVertexMap(vid), edgemodearray[n]);
2172
2173 // Set the value in the vertex edge/face matrix
2174 Svef.SetValue(0, n, VertexEdgeFaceValue);
2175 }
2176
2177 // Vertex-face coupling
2178 for (n = 0; n < nfacemodesconnected; ++n)
2179 {
2180 // Matrix value for each coefficient location
2181 VertexEdgeFaceValue = (*r_bnd)(GetVertexMap(vid), facemodearray[n]);
2182
2183 // Set the value in the vertex edge/face matrix
2184 Svef.SetValue(0, n + nedgemodesconnected, VertexEdgeFaceValue);
2185 }
2186
2187 /*
2188 * Build the edge-face transform matrix: This matrix is
2189 * constructed from the submatrices corresponding to the couping
2190 * between the edges and faces on the attached faces/edges of a
2191 * vertex.
2192 */
2193
2194 // Allocation of matrix to store edge/face-edge/face coupling
2195 DNekMatSharedPtr edgefacecoupling =
2196 MemoryManager<DNekMat>::AllocateSharedPtr(efRow, efRow, 0.0,
2197 storage);
2198 DNekMat &Sefef = (*edgefacecoupling);
2199
2200 NekDouble EdgeEdgeValue, FaceFaceValue;
2201
2202 // Edge-edge coupling (S_{ee})
2203 for (m = 0; m < nedgemodesconnected; ++m)
2204 {
2205 for (n = 0; n < nedgemodesconnected; ++n)
2206 {
2207 // Matrix value for each coefficient location
2208 EdgeEdgeValue = (*r_bnd)(edgemodearray[n], edgemodearray[m]);
2209
2210 // Set the value in the vertex edge/face matrix
2211 Sefef.SetValue(n, m, EdgeEdgeValue);
2212 }
2213 }
2214
2215 // Face-face coupling (S_{ff})
2216 for (n = 0; n < nfacemodesconnected; ++n)
2217 {
2218 for (m = 0; m < nfacemodesconnected; ++m)
2219 {
2220 // Matrix value for each coefficient location
2221 FaceFaceValue = (*r_bnd)(facemodearray[n], facemodearray[m]);
2222 // Set the value in the vertex edge/face matrix
2223 Sefef.SetValue(nedgemodesconnected + n, nedgemodesconnected + m,
2224 FaceFaceValue);
2225 }
2226 }
2227
2228 // Edge-face coupling (S_{ef} and trans(S_{ef}))
2229 for (n = 0; n < nedgemodesconnected; ++n)
2230 {
2231 for (m = 0; m < nfacemodesconnected; ++m)
2232 {
2233 // Matrix value for each coefficient location
2234 FaceFaceValue = (*r_bnd)(edgemodearray[n], facemodearray[m]);
2235
2236 // Set the value in the vertex edge/face matrix
2237 Sefef.SetValue(n, nedgemodesconnected + m, FaceFaceValue);
2238
2239 FaceFaceValue = (*r_bnd)(facemodearray[m], edgemodearray[n]);
2240
2241 // and transpose
2242 Sefef.SetValue(nedgemodesconnected + m, n, FaceFaceValue);
2243 }
2244 }
2245
2246 // Invert edge-face coupling matrix
2247 if (efRow)
2248 {
2249 Sefef.Invert();
2250
2251 // R_{v}=-S_{v,ef}inv(S_{ef,ef})
2252 Sveft = -Svef * Sefef;
2253 }
2254
2255 // Populate R with R_{ve} components
2256 for (n = 0; n < edgemodearray.size(); ++n)
2257 {
2258 R.SetValue(GetVertexMap(vid), edgemodearray[n], Sveft(0, n));
2259 }
2260
2261 // Populate R with R_{vf} components
2262 for (n = 0; n < facemodearray.size(); ++n)
2263 {
2264 R.SetValue(GetVertexMap(vid), facemodearray[n],
2265 Sveft(0, n + nedgemodesconnected));
2266 }
2267 }
2268
2269 /********************/
2270 /* edge-face matrix */
2271 /********************/
2272
2273 /*
2274 * The matrix component of \f$\mathbf{R}\f$ is given by \f[
2275 * \mathbf{R^{T}_{ef}}=-\mathbf{S}^{-1}_{ff}\mathbf{S}^{T}_{ef}\f]
2276 *
2277 * For each edge extract the submatrices from statically condensed
2278 * matrix \f$\mathbf{S}\f$ corresponding to inner products of modes
2279 * on the two attached faces within themselves as well as the
2280 * coupling matrix between the two faces
2281 * (\f$\mathbf{S}_{ff}\f$). This matrix of face coupling is then
2282 * inverted and multiplied by the submatrices of corresponding to
2283 * the coupling between the edge and attached faces
2284 * (\f$\mathbf{S}_{ef}\f$).
2285 */
2286
2287 NekDouble EdgeFaceValue, FaceFaceValue;
2288 int efCol, efRow, nedgemodes;
2289
2290 // Number of attached faces is always 2
2291 nConnectedFaces = 2;
2292
2293 // Location in the matrix
2294 MatEdgeLocation = Array<OneD, Array<OneD, unsigned int>>(nEdges);
2295 MatFaceLocation = Array<OneD, Array<OneD, unsigned int>>(nConnectedFaces);
2296
2297 // Build the edge/face transform matrix: This matrix is constructed
2298 // from the submatrices corresponding to the couping between a
2299 // specific edge and the two attached faces.
2300 for (eid = 0; eid < nEdges; ++eid)
2301 {
2302 // Row and column size of the vertex-edge/face matrix
2303 efCol = GetTraceIntNcoeffs(geom->GetEdgeFaceMap(eid, 0)) +
2304 GetTraceIntNcoeffs(geom->GetEdgeFaceMap(eid, 1));
2305 efRow = GetEdgeNcoeffs(eid) - 2;
2306
2307 // Edge-face coupling matrix
2308 DNekMatSharedPtr efedgefacecoupling =
2309 MemoryManager<DNekMat>::AllocateSharedPtr(efRow, efCol, 0.0,
2310 storage);
2311 DNekMat &Mef = (*efedgefacecoupling);
2312
2313 // Face-face coupling matrix
2314 DNekMatSharedPtr effacefacecoupling =
2315 MemoryManager<DNekMat>::AllocateSharedPtr(efCol, efCol, 0.0,
2316 storage);
2317 DNekMat &Meff = (*effacefacecoupling);
2318
2319 // Edge-face transformation matrix
2320 DNekMatSharedPtr edgefacetransformmatrix =
2321 MemoryManager<DNekMat>::AllocateSharedPtr(efRow, efCol, 0.0,
2322 storage);
2323 DNekMat &Meft = (*edgefacetransformmatrix);
2324
2325 int nfacemodesconnected =
2326 GetTraceIntNcoeffs(geom->GetEdgeFaceMap(eid, 0)) +
2327 GetTraceIntNcoeffs(geom->GetEdgeFaceMap(eid, 1));
2328 Array<OneD, unsigned int> facemodearray(nfacemodesconnected);
2329
2330 // Create array of edge modes
2331 Array<OneD, unsigned int> inedgearray = GetEdgeInverseBoundaryMap(eid);
2332 nedgemodes = GetEdgeNcoeffs(eid) - 2;
2333 Array<OneD, unsigned int> edgemodearray(nedgemodes);
2334
2335 if (nedgemodes)
2336 {
2337 Vmath::Vcopy(nedgemodes, &inedgearray[0], 1, &edgemodearray[0], 1);
2338 }
2339
2340 int offset = 0;
2341
2342 // Create array of face modes
2343 for (fid = 0; fid < nConnectedFaces; ++fid)
2344 {
2345 MatFaceLocation[fid] =
2346 GetTraceInverseBoundaryMap(geom->GetEdgeFaceMap(eid, fid));
2347 nmodes = MatFaceLocation[fid].size();
2348
2349 if (nmodes)
2350 {
2351 Vmath::Vcopy(nmodes, &MatFaceLocation[fid][0], 1,
2352 &facemodearray[offset], 1);
2353 }
2354 offset += nmodes;
2355 }
2356
2357 // Edge-face coupling
2358 for (n = 0; n < nedgemodes; ++n)
2359 {
2360 for (m = 0; m < nfacemodesconnected; ++m)
2361 {
2362 // Matrix value for each coefficient location
2363 EdgeFaceValue = (*r_bnd)(edgemodearray[n], facemodearray[m]);
2364
2365 // Set the value in the edge/face matrix
2366 Mef.SetValue(n, m, EdgeFaceValue);
2367 }
2368 }
2369
2370 // Face-face coupling
2371 for (n = 0; n < nfacemodesconnected; ++n)
2372 {
2373 for (m = 0; m < nfacemodesconnected; ++m)
2374 {
2375 // Matrix value for each coefficient location
2376 FaceFaceValue = (*r_bnd)(facemodearray[n], facemodearray[m]);
2377
2378 // Set the value in the vertex edge/face matrix
2379 Meff.SetValue(n, m, FaceFaceValue);
2380 }
2381 }
2382
2383 if (efCol)
2384 {
2385 // Invert edge-face coupling matrix
2386 Meff.Invert();
2387
2388 // trans(R_{ef})=-S_{ef}*(inv(S_{ff})
2389 Meft = -Mef * Meff;
2390 }
2391
2392 // Populate transformation matrix with Meft
2393 for (n = 0; n < Meft.GetRows(); ++n)
2394 {
2395 for (m = 0; m < Meft.GetColumns(); ++m)
2396 {
2397 R.SetValue(edgemodearray[n], facemodearray[m], Meft(n, m));
2398 }
2399 }
2400 }
2401
2402 for (i = 0; i < R.GetRows(); ++i)
2403 {
2404 R.SetValue(i, i, 1.0);
2405 }
2406
2407 if ((matrixType == StdRegions::ePreconR) ||
2408 (matrixType == StdRegions::ePreconRMass))
2409 {
2410 return transformationmatrix;
2411 }
2412 else
2413 {
2414 NEKERROR(ErrorUtil::efatal, "unkown matrix type");
2415 return NullDNekMatSharedPtr;
2416 }
2417}
2418
2419/**
2420 * \brief Build inverse and inverse transposed transformation matrix:
2421 * \f$\mathbf{R^{-1}}\f$ and \f$\mathbf{R^{-T}}\f$
2422 *
2423 * \f\mathbf{R^{-T}}=[\left[\begin{array}{ccc} \mathbf{I} &
2424 * -\mathbf{R}_{ef} & -\mathbf{R}_{ve}+\mathbf{R}_{ve}\mathbf{R}_{vf} \\
2425 * 0 & \mathbf{I} & \mathbf{R}_{ef} \\
2426 * 0 & 0 & \mathbf{I}} \end{array}\right]\f]
2427 */
2428DNekMatSharedPtr Expansion3D::v_BuildInverseTransformationMatrix(
2429 const DNekScalMatSharedPtr &transformationmatrix)
2430{
2431 int i, j, n, eid = 0, fid = 0;
2432 int nCoeffs = NumBndryCoeffs();
2433 NekDouble MatrixValue;
2434 DNekScalMat &R = (*transformationmatrix);
2435
2436 // Define storage for vertex transpose matrix and zero all entries
2437 MatrixStorage storage = eFULL;
2438
2439 // Inverse transformation matrix
2440 DNekMatSharedPtr inversetransformationmatrix =
2441 MemoryManager<DNekMat>::AllocateSharedPtr(nCoeffs, nCoeffs, 0.0,
2442 storage);
2443 DNekMat &InvR = (*inversetransformationmatrix);
2444
2445 int nVerts = GetNverts();
2446 int nEdges = GetNedges();
2447 int nFaces = GetNtraces();
2448
2449 int nedgemodes = 0;
2450 int nfacemodes = 0;
2451 int nedgemodestotal = 0;
2452 int nfacemodestotal = 0;
2453
2454 for (eid = 0; eid < nEdges; ++eid)
2455 {
2456 nedgemodes = GetEdgeNcoeffs(eid) - 2;
2457 nedgemodestotal += nedgemodes;
2458 }
2459
2460 for (fid = 0; fid < nFaces; ++fid)
2461 {
2462 nfacemodes = GetTraceIntNcoeffs(fid);
2463 nfacemodestotal += nfacemodes;
2464 }
2465
2466 Array<OneD, unsigned int> edgemodearray(nedgemodestotal);
2467 Array<OneD, unsigned int> facemodearray(nfacemodestotal);
2468
2469 int offset = 0;
2470
2471 // Create array of edge modes
2472 for (eid = 0; eid < nEdges; ++eid)
2473 {
2474 Array<OneD, unsigned int> edgearray = GetEdgeInverseBoundaryMap(eid);
2475 nedgemodes = GetEdgeNcoeffs(eid) - 2;
2476
2477 // Only copy if there are edge modes
2478 if (nedgemodes)
2479 {
2480 Vmath::Vcopy(nedgemodes, &edgearray[0], 1, &edgemodearray[offset],
2481 1);
2482 }
2483
2484 offset += nedgemodes;
2485 }
2486
2487 offset = 0;
2488
2489 // Create array of face modes
2490 for (fid = 0; fid < nFaces; ++fid)
2491 {
2492 Array<OneD, unsigned int> facearray = GetTraceInverseBoundaryMap(fid);
2493 nfacemodes = GetTraceIntNcoeffs(fid);
2494
2495 // Only copy if there are face modes
2496 if (nfacemodes)
2497 {
2498 Vmath::Vcopy(nfacemodes, &facearray[0], 1, &facemodearray[offset],
2499 1);
2500 }
2501
2502 offset += nfacemodes;
2503 }
2504
2505 // Vertex-edge/face
2506 for (i = 0; i < nVerts; ++i)
2507 {
2508 for (j = 0; j < nedgemodestotal; ++j)
2509 {
2510 InvR.SetValue(GetVertexMap(i), edgemodearray[j],
2511 -R(GetVertexMap(i), edgemodearray[j]));
2512 }
2513 for (j = 0; j < nfacemodestotal; ++j)
2514 {
2515 InvR.SetValue(GetVertexMap(i), facemodearray[j],
2516 -R(GetVertexMap(i), facemodearray[j]));
2517 for (n = 0; n < nedgemodestotal; ++n)
2518 {
2519 MatrixValue = InvR.GetValue(GetVertexMap(i), facemodearray[j]) +
2520 R(GetVertexMap(i), edgemodearray[n]) *
2521 R(edgemodearray[n], facemodearray[j]);
2522 InvR.SetValue(GetVertexMap(i), facemodearray[j], MatrixValue);
2523 }
2524 }
2525 }
2526
2527 // Edge-face contributions
2528 for (i = 0; i < nedgemodestotal; ++i)
2529 {
2530 for (j = 0; j < nfacemodestotal; ++j)
2531 {
2532 InvR.SetValue(edgemodearray[i], facemodearray[j],
2533 -R(edgemodearray[i], facemodearray[j]));
2534 }
2535 }
2536
2537 for (i = 0; i < nCoeffs; ++i)
2538 {
2539 InvR.SetValue(i, i, 1.0);
2540 }
2541
2542 return inversetransformationmatrix;
2543}
2544
2545Array<OneD, unsigned int> Expansion3D::GetEdgeInverseBoundaryMap(int eid)
2546{
2547 int n, j;
2548 int nEdgeCoeffs;
2549 int nBndCoeffs = NumBndryCoeffs();
2550
2551 Array<OneD, unsigned int> bmap(nBndCoeffs);
2552 GetBoundaryMap(bmap);
2553
2554 // Map from full system to statically condensed system (i.e reverse
2555 // GetBoundaryMap)
2556 map<int, int> invmap;
2557 for (j = 0; j < nBndCoeffs; ++j)
2558 {
2559 invmap[bmap[j]] = j;
2560 }
2561
2562 // Number of interior edge coefficients
2563 nEdgeCoeffs = GetEdgeNcoeffs(eid) - 2;
2564
2565 const SpatialDomains::Geometry3D *geom = GetGeom3D();
2566
2567 Array<OneD, unsigned int> edgemaparray(nEdgeCoeffs);
2568 StdRegions::Orientation eOrient = geom->GetEorient(eid);
2569 Array<OneD, unsigned int> maparray = Array<OneD, unsigned int>(nEdgeCoeffs);
2570 Array<OneD, int> signarray = Array<OneD, int>(nEdgeCoeffs, 1);
2571
2572 // maparray is the location of the edge within the matrix
2573 GetEdgeInteriorToElementMap(eid, maparray, signarray, eOrient);
2574
2575 for (n = 0; n < nEdgeCoeffs; ++n)
2576 {
2577 edgemaparray[n] = invmap[maparray[n]];
2578 }
2579
2580 return edgemaparray;
2581}
2582
2583Array<OneD, unsigned int> Expansion3D::GetTraceInverseBoundaryMap(
2584 int fid, StdRegions::Orientation faceOrient, int P1, int P2)
2585{
2586 int n, j;
2587 int nFaceCoeffs;
2588
2589 int nBndCoeffs = NumBndryCoeffs();
2590
2591 Array<OneD, unsigned int> bmap(nBndCoeffs);
2592 GetBoundaryMap(bmap);
2593
2594 // Map from full system to statically condensed system (i.e reverse
2595 // GetBoundaryMap)
2596 map<int, int> reversemap;
2597 for (j = 0; j < bmap.size(); ++j)
2598 {
2599 reversemap[bmap[j]] = j;
2600 }
2601
2602 // Number of interior face coefficients
2603 nFaceCoeffs = GetTraceIntNcoeffs(fid);
2604
2605 StdRegions::Orientation fOrient;
2606 Array<OneD, unsigned int> maparray = Array<OneD, unsigned int>(nFaceCoeffs);
2607 Array<OneD, int> signarray = Array<OneD, int>(nFaceCoeffs, 1);
2608
2609 if (faceOrient == StdRegions::eNoOrientation)
2610 {
2611 fOrient = GetTraceOrient(fid);
2612 }
2613 else
2614 {
2615 fOrient = faceOrient;
2616 }
2617
2618 // maparray is the location of the face within the matrix
2619 GetTraceInteriorToElementMap(fid, maparray, signarray, fOrient);
2620
2621 Array<OneD, unsigned int> facemaparray;
2622 int locP1, locP2;
2623 GetTraceNumModes(fid, locP1, locP2, fOrient);
2624
2625 if (P1 == -1)
2626 {
2627 P1 = locP1;
2628 }
2629 else
2630 {
2631 ASSERTL1(P1 <= locP1, "Expect value of passed P1 to "
2632 "be lower or equal to face num modes");
2633 }
2634
2635 if (P2 == -1)
2636 {
2637 P2 = locP2;
2638 }
2639 else
2640 {
2641 ASSERTL1(P2 <= locP2, "Expect value of passed P2 to "
2642 "be lower or equal to face num modes");
2643 }
2644
2645 switch (GetGeom3D()->GetFace(fid)->GetShapeType())
2646 {
2647 case LibUtilities::eTriangle:
2648 {
2649 if (((P1 - 3) > 0) && ((P2 - 3) > 0))
2650 {
2651 facemaparray = Array<OneD, unsigned int>(
2652 LibUtilities::StdTriData::getNumberOfCoefficients(P1 - 3,
2653 P2 - 3));
2654 int cnt = 0;
2655 int cnt1 = 0;
2656 for (n = 0; n < P1 - 3; ++n)
2657 {
2658 for (int m = 0; m < P2 - 3 - n; ++m, ++cnt)
2659 {
2660 facemaparray[cnt] = reversemap[maparray[cnt1 + m]];
2661 }
2662 cnt1 += locP2 - 3 - n;
2663 }
2664 }
2665 }
2666 break;
2667 case LibUtilities::eQuadrilateral:
2668 {
2669 if (((P1 - 2) > 0) && ((P2 - 2) > 0))
2670 {
2671 facemaparray = Array<OneD, unsigned int>(
2672 LibUtilities::StdQuadData::getNumberOfCoefficients(P1 - 2,
2673 P2 - 2));
2674 int cnt = 0;
2675 int cnt1 = 0;
2676 for (n = 0; n < P2 - 2; ++n)
2677 {
2678 for (int m = 0; m < P1 - 2; ++m, ++cnt)
2679 {
2680 facemaparray[cnt] = reversemap[maparray[cnt1 + m]];
2681 }
2682 cnt1 += locP1 - 2;
2683 }
2684 }
2685 }
2686 break;
2687 default:
2688 {
2689 ASSERTL0(false, "Invalid shape type.");
2690 }
2691 break;
2692 }
2693
2694 return facemaparray;
2695}
2696
2697void Expansion3D::GetInverseBoundaryMaps(
2698 Array<OneD, unsigned int> &vmap,
2699 Array<OneD, Array<OneD, unsigned int>> &emap,
2700 Array<OneD, Array<OneD, unsigned int>> &fmap)
2701{
2702 int n, j;
2703 int nEdgeCoeffs;
2704 int nFaceCoeffs;
2705
2706 int nBndCoeffs = NumBndryCoeffs();
2707
2708 Array<OneD, unsigned int> bmap(nBndCoeffs);
2709 GetBoundaryMap(bmap);
2710
2711 // Map from full system to statically condensed system (i.e reverse
2712 // GetBoundaryMap)
2713 map<int, int> reversemap;
2714 for (j = 0; j < bmap.size(); ++j)
2715 {
2716 reversemap[bmap[j]] = j;
2717 }
2718
2719 int nverts = GetNverts();
2720 vmap = Array<OneD, unsigned int>(nverts);
2721 for (n = 0; n < nverts; ++n)
2722 {
2723 int id = GetVertexMap(n);
2724 vmap[n] = reversemap[id]; // not sure what should be true here.
2725 }
2726
2727 int nedges = GetNedges();
2728 emap = Array<OneD, Array<OneD, unsigned int>>(nedges);
2729
2730 for (int eid = 0; eid < nedges; ++eid)
2731 {
2732 // Number of interior edge coefficients
2733 nEdgeCoeffs = GetEdgeNcoeffs(eid) - 2;
2734
2735 Array<OneD, unsigned int> edgemaparray(nEdgeCoeffs);
2736 Array<OneD, unsigned int> maparray =
2737 Array<OneD, unsigned int>(nEdgeCoeffs);
2738 Array<OneD, int> signarray = Array<OneD, int>(nEdgeCoeffs, 1);
2739
2740 // maparray is the location of the edge within the matrix
2741 GetEdgeInteriorToElementMap(eid, maparray, signarray,
2742 StdRegions::eForwards);
2743
2744 for (n = 0; n < nEdgeCoeffs; ++n)
2745 {
2746 edgemaparray[n] = reversemap[maparray[n]];
2747 }
2748 emap[eid] = edgemaparray;
2749 }
2750
2751 int nfaces = GetNtraces();
2752 fmap = Array<OneD, Array<OneD, unsigned int>>(nfaces);
2753
2754 for (int fid = 0; fid < nfaces; ++fid)
2755 {
2756 // Number of interior face coefficients
2757 nFaceCoeffs = GetTraceIntNcoeffs(fid);
2758
2759 Array<OneD, unsigned int> facemaparray(nFaceCoeffs);
2760 Array<OneD, unsigned int> maparray =
2761 Array<OneD, unsigned int>(nFaceCoeffs);
2762 Array<OneD, int> signarray = Array<OneD, int>(nFaceCoeffs, 1);
2763
2764 // maparray is the location of the face within the matrix
2765 GetTraceInteriorToElementMap(fid, maparray, signarray,
2766 StdRegions::eDir1FwdDir1_Dir2FwdDir2);
2767
2768 for (n = 0; n < nFaceCoeffs; ++n)
2769 {
2770 facemaparray[n] = reversemap[maparray[n]];
2771 }
2772
2773 fmap[fid] = facemaparray;
2774 }
2775}
2776
2781
2782/**
2783 * @brief Extract the physical values along face \a face from \a
2784 * inarray into \a outarray following the local face orientation
2785 * and point distribution defined by defined in \a FaceExp.
2786 */
2787void Expansion3D::v_GetTracePhysVals(
2788 const int face, const StdRegions::StdExpansionSharedPtr &FaceExp,
2789 const Array<OneD, const NekDouble> &inarray,
2790 Array<OneD, NekDouble> &outarray, StdRegions::Orientation orient)
2791{
2792 if (orient == StdRegions::eNoOrientation)
2793 {
2794 orient = GetTraceOrient(face);
2795 }
2796
2797 int nq0 = FaceExp->GetNumPoints(0);
2798 int nq1 = FaceExp->GetNumPoints(1);
2799
2800 int nfacepts = GetTraceNumPoints(face);
2801 int dir0 = GetGeom3D()->GetDir(face, 0);
2802 int dir1 = GetGeom3D()->GetDir(face, 1);
2803
2804 Array<OneD, NekDouble> o_tmp(nfacepts);
2805 Array<OneD, NekDouble> o_tmp2(FaceExp->GetTotPoints());
2806 Array<OneD, int> faceids;
2807
2808 // Get local face pts and put into o_tmp
2809 GetTracePhysMap(face, faceids);
2810 // The static cast is necessary because faceids should be
2811 // Array<OneD, size_t> faceids ... or at leaste the same type as
2812 // faceids.size() ...
2813 Vmath::Gathr(static_cast<int>(faceids.size()), inarray, faceids, o_tmp);
2814
2815 int to_id0, to_id1;
2816
2817 if (orient < StdRegions::eDir1FwdDir2_Dir2FwdDir1)
2818 {
2819 to_id0 = 0;
2820 to_id1 = 1;
2821 }
2822 else // transpose points key evaluation
2823 {
2824 to_id0 = 1;
2825 to_id1 = 0;
2826 }
2827
2828 // interpolate to points distrbution given in FaceExp
2829 LibUtilities::Interp2D(
2830 m_base[dir0]->GetPointsKey(), m_base[dir1]->GetPointsKey(),
2831 o_tmp.data(), FaceExp->GetBasis(to_id0)->GetPointsKey(),
2832 FaceExp->GetBasis(to_id1)->GetPointsKey(), o_tmp2.data());
2833
2834 // Reshuffule points as required and put into outarray.
2835 v_ReOrientTracePhysMap(orient, faceids, nq0, nq1);
2836 Vmath::Scatr(nq0 * nq1, o_tmp2, faceids, outarray);
2837}
2838
2839void Expansion3D::v_GenTraceExp(const int traceid, ExpansionSharedPtr &exp)
2840{
2841 SpatialDomains::Geometry *faceGeom = m_geom->GetFace(traceid);
2842 if (faceGeom->GetNumVerts() == 3)
2843 {
2844 exp = MemoryManager<LocalRegions::TriExp>::AllocateSharedPtr(
2845 GetTraceBasisKey(traceid, 0), GetTraceBasisKey(traceid, 1),
2846 m_geom->GetFace(traceid));
2847 }
2848 else
2849 {
2850 exp = MemoryManager<LocalRegions::QuadExp>::AllocateSharedPtr(
2851 GetTraceBasisKey(traceid, 0), GetTraceBasisKey(traceid, 1),
2852 m_geom->GetFace(traceid));
2853 }
2854}
2855
2856void Expansion3D::v_ReOrientTracePhysMap(const StdRegions::Orientation orient,
2857 Array<OneD, int> &idmap, const int nq0,
2858 const int nq1)
2859{
2860
2861 if (idmap.size() != nq0 * nq1)
2862 {
2863 idmap = Array<OneD, int>(nq0 * nq1);
2864 }
2865
2866 if (GetNverts() == 3) // Tri face
2867 {
2868 switch (orient)
2869 {
2870 case StdRegions::eDir1FwdDir1_Dir2FwdDir2:
2871 // eseentially straight copy
2872 for (int i = 0; i < nq0 * nq1; ++i)
2873 {
2874 idmap[i] = i;
2875 }
2876 break;
2877 case StdRegions::eDir1BwdDir1_Dir2FwdDir2:
2878 // reverse
2879 for (int j = 0; j < nq1; ++j)
2880 {
2881 for (int i = 0; i < nq0; ++i)
2882 {
2883 idmap[j * nq0 + i] = nq0 - 1 - i + j * nq0;
2884 }
2885 }
2886 break;
2887 default:
2888 ASSERTL0(false,
2889 "Case not supposed to be used in this function");
2890 }
2891 }
2892 else
2893 {
2894 switch (orient)
2895 {
2896 case StdRegions::eDir1FwdDir1_Dir2FwdDir2:
2897 // eseentially straight copy
2898 for (int i = 0; i < nq0 * nq1; ++i)
2899 {
2900 idmap[i] = i;
2901 }
2902 break;
2903 case StdRegions::eDir1BwdDir1_Dir2FwdDir2:
2904 {
2905 // Direction A negative and B positive
2906 for (int j = 0; j < nq1; j++)
2907 {
2908 for (int i = 0; i < nq0; ++i)
2909 {
2910 idmap[j * nq0 + i] = nq0 - 1 - i + j * nq0;
2911 }
2912 }
2913 }
2914 break;
2915 case StdRegions::eDir1FwdDir1_Dir2BwdDir2:
2916 {
2917 // Direction A positive and B negative
2918 for (int j = 0; j < nq1; j++)
2919 {
2920 for (int i = 0; i < nq0; ++i)
2921 {
2922 idmap[j * nq0 + i] = nq0 * (nq1 - 1 - j) + i;
2923 }
2924 }
2925 }
2926 break;
2927 case StdRegions::eDir1BwdDir1_Dir2BwdDir2:
2928 {
2929 // Direction A negative and B negative
2930 for (int j = 0; j < nq1; j++)
2931 {
2932 for (int i = 0; i < nq0; ++i)
2933 {
2934 idmap[j * nq0 + i] = nq0 * nq1 - 1 - j * nq0 - i;
2935 }
2936 }
2937 }
2938 break;
2939 case StdRegions::eDir1FwdDir2_Dir2FwdDir1:
2940 {
2941 // Transposed, Direction A and B positive
2942 for (int i = 0; i < nq0; ++i)
2943 {
2944 for (int j = 0; j < nq1; ++j)
2945 {
2946 idmap[i * nq1 + j] = i + j * nq0;
2947 }
2948 }
2949 }
2950 break;
2951 case StdRegions::eDir1FwdDir2_Dir2BwdDir1:
2952 {
2953 // Transposed, Direction A positive and B negative
2954 for (int i = 0; i < nq0; ++i)
2955 {
2956 for (int j = 0; j < nq1; ++j)
2957 {
2958 idmap[i * nq1 + j] = nq0 - 1 - i + j * nq0;
2959 }
2960 }
2961 }
2962 break;
2963 case StdRegions::eDir1BwdDir2_Dir2FwdDir1:
2964 {
2965 // Transposed, Direction A negative and B positive
2966 for (int i = 0; i < nq0; ++i)
2967 {
2968 for (int j = 0; j < nq1; ++j)
2969 {
2970 idmap[i * nq1 + j] = i + nq0 * (nq1 - 1) - j * nq0;
2971 }
2972 }
2973 }
2974 break;
2975 case StdRegions::eDir1BwdDir2_Dir2BwdDir1:
2976 {
2977 // Transposed, Direction A and B negative
2978 for (int i = 0; i < nq0; ++i)
2979 {
2980 for (int j = 0; j < nq1; ++j)
2981 {
2982 idmap[i * nq1 + j] = nq0 * nq1 - 1 - i - j * nq0;
2983 }
2984 }
2985 }
2986 break;
2987 default:
2988 ASSERTL0(false, "Unknow orientation");
2989 break;
2990 }
2991 }
2992}
2993
2994void Expansion3D::v_NormVectorIProductWRTBase(
2995 const Array<OneD, const Array<OneD, NekDouble>> &Fvec,
2996 Array<OneD, NekDouble> &outarray)
2997{
2998 NormVectorIProductWRTBase(Fvec[0], Fvec[1], Fvec[2], outarray);
2999}
3000
3001// Compute edgenormal \cdot vector
3002Array<OneD, NekDouble> Expansion3D::GetnFacecdotMF(
3003 const int dir, const int face, ExpansionSharedPtr &FaceExp_f,
3004 const Array<OneD, const Array<OneD, NekDouble>> &normals,
3005 const StdRegions::VarCoeffMap &varcoeffs)
3006{
3007 int nquad_f = FaceExp_f->GetNumPoints(0) * FaceExp_f->GetNumPoints(1);
3008 int coordim = GetCoordim();
3009
3010 int nquad0 = m_base[0]->GetNumPoints();
3011 int nquad1 = m_base[1]->GetNumPoints();
3012 int nquad2 = m_base[2]->GetNumPoints();
3013 int nqtot = nquad0 * nquad1 * nquad2;
3014
3015 StdRegions::VarCoeffType MMFCoeffs[15] = {
3016 StdRegions::eVarCoeffMF1x, StdRegions::eVarCoeffMF1y,
3017 StdRegions::eVarCoeffMF1z, StdRegions::eVarCoeffMF1Div,
3018 StdRegions::eVarCoeffMF1Mag, StdRegions::eVarCoeffMF2x,
3019 StdRegions::eVarCoeffMF2y, StdRegions::eVarCoeffMF2z,
3020 StdRegions::eVarCoeffMF2Div, StdRegions::eVarCoeffMF2Mag,
3021 StdRegions::eVarCoeffMF3x, StdRegions::eVarCoeffMF3y,
3022 StdRegions::eVarCoeffMF3z, StdRegions::eVarCoeffMF3Div,
3023 StdRegions::eVarCoeffMF3Mag};
3024
3025 StdRegions::VarCoeffMap::const_iterator MFdir;
3026
3027 Array<OneD, NekDouble> nFacecdotMF(nqtot, 0.0);
3028 Array<OneD, NekDouble> tmp(nqtot);
3029 Array<OneD, NekDouble> tmp_f(nquad_f);
3030 for (int k = 0; k < coordim; k++)
3031 {
3032 MFdir = varcoeffs.find(MMFCoeffs[dir * 5 + k]);
3033 tmp = MFdir->second.GetValue();
3034
3035 GetPhysFaceVarCoeffsFromElement(face, FaceExp_f, tmp, tmp_f);
3036
3037 Vmath::Vvtvp(nquad_f, &tmp_f[0], 1, &normals[k][0], 1, &nFacecdotMF[0],
3038 1, &nFacecdotMF[0], 1);
3039 }
3040
3041 return nFacecdotMF;
3042}
3043
3044void Expansion3D::v_TraceNormLen(const int traceid, NekDouble &h, NekDouble &p)
3045{
3046 SpatialDomains::Geometry *geom = GetGeom();
3047
3048 int nverts = geom->GetFace(traceid)->GetNumVerts();
3049
3050 SpatialDomains::PointGeom tn1, tn2, normal;
3051 tn1.Sub(*(geom->GetFace(traceid)->GetVertex(1)),
3052 *(geom->GetFace(traceid)->GetVertex(0)));
3053
3054 tn2.Sub(*(geom->GetFace(traceid)->GetVertex(nverts - 1)),
3055 *(geom->GetFace(traceid)->GetVertex(0)));
3056
3057 normal.Mult(tn1, tn2);
3058
3059 // normalise normal
3060 NekDouble mag = normal.dot(normal);
3061 mag = 1.0 / sqrt(mag);
3062 normal.UpdatePosition(normal.x() * mag, normal.y() * mag, normal.z() * mag);
3063
3064 SpatialDomains::PointGeom Dx;
3065 h = 0.0;
3066 p = 0.0;
3067 for (int i = 0; i < nverts; ++i)
3068 {
3069 // vertices on edges
3070 int edgid = geom->GetEdgeNormalToFaceVert(traceid, i);
3071
3072 // vector along noramal edge to each vertex
3073 Dx.Sub(*(geom->GetEdge(edgid)->GetVertex(0)),
3074 *(geom->GetEdge(edgid)->GetVertex(1)));
3075
3076 // calculate perpendicular distance of normal length
3077 // from first vertex
3078 h += fabs(normal.dot(Dx));
3079 }
3080
3081 h /= static_cast<NekDouble>(nverts);
3082
3083 // find normal basis direction
3084 int dir0 = geom->GetDir(traceid, 0);
3085 int dir1 = geom->GetDir(traceid, 1);
3086 int dirn;
3087 for (dirn = 0; dirn < 3; ++dirn)
3088 {
3089 if ((dirn != dir0) && (dirn != dir1))
3090 {
3091 break;
3092 }
3093 }
3094 p = (NekDouble)(GetBasisNumModes(dirn) - 1);
3095}
3096} // namespace Nektar::LocalRegions
ASSERTL0
#define ASSERTL0(condition, msg)
Definition ErrorUtil.hpp:216
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition ErrorUtil.hpp:210
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::Expansion3D::v_BuildTransformationMatrix
DNekMatSharedPtr v_BuildTransformationMatrix(const DNekScalMatSharedPtr &r_bnd, const StdRegions::MatrixType matrixType) override
Definition Expansion3D.cpp:2046
Nektar::LocalRegions::Expansion3D::v_AddRobinMassMatrix
void v_AddRobinMassMatrix(const int face, const Array< OneD, const NekDouble > &primCoeffs, DNekMatSharedPtr &inoutmat) override
Definition Expansion3D.cpp:1862
Nektar::LocalRegions::Expansion3D::SetTraceToGeomOrientation
void SetTraceToGeomOrientation(Array< OneD, NekDouble > &inout)
Align trace orientation with the geometry orientation.
Definition Expansion3D.cpp:393
Nektar::LocalRegions::Expansion3D::v_DGDeriv
void v_DGDeriv(const int dir, const Array< OneD, const NekDouble > &incoeffs, Array< OneD, ExpansionSharedPtr > &FaceExp, Array< OneD, Array< OneD, NekDouble > > &faceCoeffs, Array< OneD, NekDouble > &out_d) override
Evaluate coefficients of weak deriviative in the direction dir given the input coefficicents incoeffs...
Definition Expansion3D.cpp:1833
Nektar::LocalRegions::Expansion3D::v_GetTraceOrient
StdRegions::Orientation v_GetTraceOrient(int face) override
Definition Expansion3D.cpp:2777
Nektar::LocalRegions::Expansion3D::AddNormTraceInt
void AddNormTraceInt(const int dir, Array< OneD, ExpansionSharedPtr > &FaceExp, Array< OneD, Array< OneD, NekDouble > > &faceCoeffs, Array< OneD, NekDouble > &outarray)
Definition Expansion3D.cpp:294
Nektar::LocalRegions::Expansion3D::CreateMatrix
DNekScalMatSharedPtr CreateMatrix(const MatrixKey &mkey)
Definition Expansion3D.cpp:407
Nektar::LocalRegions::Expansion3D::GetTraceInverseBoundaryMap
Array< OneD, unsigned int > GetTraceInverseBoundaryMap(int fid, StdRegions::Orientation faceOrient=StdRegions::eNoOrientation, int P1=-1, int P2=-1)
Definition Expansion3D.cpp:2583
Nektar::LocalRegions::Expansion3D::m_requireNeg
std::vector< bool > m_requireNeg
Definition Expansion3D.h:168
Nektar::LocalRegions::Expansion3D::GetInverseBoundaryMaps
void GetInverseBoundaryMaps(Array< OneD, unsigned int > &vmap, Array< OneD, Array< OneD, unsigned int > > &emap, Array< OneD, Array< OneD, unsigned int > > &fmap)
Definition Expansion3D.cpp:2697
Nektar::LocalRegions::Expansion3D::GetnFacecdotMF
Array< OneD, NekDouble > GetnFacecdotMF(const int dir, const int face, ExpansionSharedPtr &FaceExp_f, const Array< OneD, const Array< OneD, NekDouble > > &normals, const StdRegions::VarCoeffMap &varcoeffs)
Definition Expansion3D.cpp:3002
Nektar::LocalRegions::Expansion3D::v_BuildInverseTransformationMatrix
DNekMatSharedPtr v_BuildInverseTransformationMatrix(const DNekScalMatSharedPtr &transformationmatrix) override
Build inverse and inverse transposed transformation matrix: and .
Definition Expansion3D.cpp:2428
Nektar::LocalRegions::Expansion3D::AddFaceBoundaryInt
void AddFaceBoundaryInt(const int face, ExpansionSharedPtr &FaceExp, Array< OneD, NekDouble > &facePhys, Array< OneD, NekDouble > &outarray, const StdRegions::VarCoeffMap &varcoeffs=StdRegions::NullVarCoeffMap)
Definition Expansion3D.cpp:322
Nektar::LocalRegions::Expansion3D::GetPhysFaceVarCoeffsFromElement
void GetPhysFaceVarCoeffsFromElement(const int face, ExpansionSharedPtr &FaceExp, const Array< OneD, const NekDouble > &varcoeff, Array< OneD, NekDouble > &outarray)
Definition Expansion3D.cpp:202
Nektar::LocalRegions::Expansion3D::v_TraceNormLen
void v_TraceNormLen(const int traceid, NekDouble &h, NekDouble &p) override
Definition Expansion3D.cpp:3044
Nektar::LocalRegions::Expansion3D::GetGeom3D
SpatialDomains::Geometry3D * GetGeom3D() const
Definition Expansion3D.h:176
Nektar::LocalRegions::Expansion3D::v_ReOrientTracePhysMap
void v_ReOrientTracePhysMap(const StdRegions::Orientation orient, Array< OneD, int > &idmap, const int nq0, const int nq1) override
Definition Expansion3D.cpp:2856
Nektar::LocalRegions::Expansion3D::v_NormVectorIProductWRTBase
void v_NormVectorIProductWRTBase(const Array< OneD, const Array< OneD, NekDouble > > &Fvec, Array< OneD, NekDouble > &outarray) override
Definition Expansion3D.cpp:2994
Nektar::LocalRegions::Expansion3D::GetEdgeInverseBoundaryMap
Array< OneD, unsigned int > GetEdgeInverseBoundaryMap(int eid)
Definition Expansion3D.cpp:2545
Nektar::LocalRegions::Expansion3D::v_GenMatrix
DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
Definition Expansion3D.cpp:936
Nektar::LocalRegions::Expansion3D::v_GenTraceExp
void v_GenTraceExp(const int traceid, ExpansionSharedPtr &exp) override
Definition Expansion3D.cpp:2839
Nektar::LocalRegions::Expansion3D::v_AddFaceNormBoundaryInt
void v_AddFaceNormBoundaryInt(const int face, const ExpansionSharedPtr &FaceExp, const Array< OneD, const NekDouble > &Fn, Array< OneD, NekDouble > &outarray) override
Definition Expansion3D.cpp:1705
Nektar::LocalRegions::Expansion3D::AddHDGHelmholtzFaceTerms
void AddHDGHelmholtzFaceTerms(const NekDouble tau, const int edge, Array< OneD, NekDouble > &facePhys, const StdRegions::VarCoeffMap &dirForcing, Array< OneD, NekDouble > &outarray)
Definition Expansion3D.cpp:50
Nektar::LocalRegions::Expansion3D::v_BuildVertexMatrix
DNekMatSharedPtr v_BuildVertexMatrix(const DNekScalMatSharedPtr &r_bnd) override
Definition Expansion3D.cpp:2016
Nektar::LocalRegions::Expansion3D::v_GetTracePhysVals
void v_GetTracePhysVals(const int face, const StdRegions::StdExpansionSharedPtr &FaceExp, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, StdRegions::Orientation orient) override
Extract the physical values along face face from inarray into outarray following the local face orien...
Definition Expansion3D.cpp:2787
Nektar::LocalRegions::Expansion3D::SetFaceToGeomOrientation
void SetFaceToGeomOrientation(const int face, Array< OneD, NekDouble > &inout)
Align face orientation with the geometry orientation.
Definition Expansion3D.cpp:348
Nektar::LocalRegions::Expansion::GetTracePhysMap
void GetTracePhysMap(const int edge, Array< OneD, int > &outarray)
Definition Expansion.h:209
Nektar::LocalRegions::Expansion::DropLocMatrix
void DropLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition Expansion.cpp:90
Nektar::LocalRegions::Expansion::GetGeom
SpatialDomains::Geometry * GetGeom() const
Definition Expansion.cpp:167
Nektar::LocalRegions::Expansion::BuildVertexMatrix
DNekMatSharedPtr BuildVertexMatrix(const DNekScalMatSharedPtr &r_bnd)
Definition Expansion.cpp:101
Nektar::LocalRegions::Expansion::m_geom
SpatialDomains::Geometry * m_geom
Definition Expansion.h:279
Nektar::LocalRegions::Expansion::GetMFMag
Array< OneD, NekDouble > GetMFMag(const int dir, const StdRegions::VarCoeffMap &varcoeffs)
Definition Expansion.cpp:794
Nektar::LocalRegions::Expansion::GetTracePhysVals
void GetTracePhysVals(const int trace, const StdRegions::StdExpansionSharedPtr &TraceExp, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, StdRegions::Orientation orient=StdRegions::eNoOrientation)
Definition Expansion.h:200
Nektar::LocalRegions::Expansion::m_traceExp
std::map< int, ExpansionWeakPtr > m_traceExp
Definition Expansion.h:278
Nektar::LocalRegions::Expansion::m_metricinfo
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
Definition Expansion.h:280
Nektar::LocalRegions::Expansion::GetTraceExp
ExpansionSharedPtr GetTraceExp(const int traceid)
Definition Expansion.h:420
Nektar::LocalRegions::Expansion::GetTraceOrient
StdRegions::Orientation GetTraceOrient(int trace)
Definition Expansion.h:168
Nektar::LocalRegions::Expansion::GetIndexMap
IndexMapValuesSharedPtr GetIndexMap(const IndexMapKey &ikey)
Definition Expansion.h:146
Nektar::LocalRegions::Expansion::GetLocMatrix
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition Expansion.cpp:84
Nektar::LocalRegions::Expansion::GetMFDiv
Array< OneD, NekDouble > GetMFDiv(const int dir, const StdRegions::VarCoeffMap &varcoeffs)
Definition Expansion.cpp:771
Nektar::LocalRegions::Expansion::TraceNormLen
void TraceNormLen(const int traceid, NekDouble &h, NekDouble &p)
Definition Expansion.h:250
Nektar::LocalRegions::Expansion::GetTraceNormal
const NormalVector & GetTraceNormal(const int id)
Definition Expansion.cpp:251
Nektar::LocalRegions::Expansion::GetMF
Array< OneD, NekDouble > GetMF(const int dir, const int shapedim, const StdRegions::VarCoeffMap &varcoeffs)
Definition Expansion.cpp:718
Nektar::LocalRegions::Expansion::BuildTransformationMatrix
DNekMatSharedPtr BuildTransformationMatrix(const DNekScalMatSharedPtr &r_bnd, const StdRegions::MatrixType matrixType)
Definition Expansion.cpp:95
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition NekMemoryManager.hpp:166
Nektar::NekPoint::x
boost::call_traits< DataType >::const_reference x() const
Definition NekPoint.hpp:160
Nektar::NekPoint::z
boost::call_traits< DataType >::const_reference z() const
Definition NekPoint.hpp:172
Nektar::NekPoint::y
boost::call_traits< DataType >::const_reference y() const
Definition NekPoint.hpp:166
Nektar::SpatialDomains::Geometry3D
3D geometry information
Definition Geometry3D.h:50
Nektar::SpatialDomains::Geometry
Base class for shape geometry information.
Definition Geometry.h:74
Nektar::SpatialDomains::Geometry::GetGlobalID
int GetGlobalID(void) const
Get the ID of this object.
Definition Geometry.h:322
Nektar::SpatialDomains::Geometry::GetVertex
PointGeom * GetVertex(int i) const
Returns vertex i of this object.
Definition Geometry.h:361
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:661
Nektar::SpatialDomains::Geometry::GetVertexEdgeMap
int GetVertexEdgeMap(int i, int j) const
Returns the standard element edge IDs that are connected to a given vertex.
Definition Geometry.h:587
Nektar::SpatialDomains::Geometry::GetNumVerts
int GetNumVerts() const
Get the number of vertices of this object.
Definition Geometry.h:403
Nektar::SpatialDomains::Geometry::GetVertexFaceMap
int GetVertexFaceMap(int i, int j) const
Returns the standard element face IDs that are connected to a given vertex.
Definition Geometry.h:609
Nektar::SpatialDomains::Geometry::GetEdge
Geometry1D * GetEdge(int i) const
Returns edge i of this object.
Definition Geometry.h:369
Nektar::SpatialDomains::Geometry::GetFace
Geometry2D * GetFace(int i) const
Returns face i of this object.
Definition Geometry.h:377
Nektar::SpatialDomains::Geometry::GetEdgeFaceMap
int GetEdgeFaceMap(int i, int j) const
Returns the standard element edge IDs that are connected to a given face.
Definition Geometry.h:630
Nektar::SpatialDomains::Geometry::GetEdgeNormalToFaceVert
int GetEdgeNormalToFaceVert(int i, int j) const
Returns the standard lement edge IDs that are normal to a given face vertex.
Definition Geometry.h:652
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:386
Nektar::SpatialDomains::Geometry::GetForient
StdRegions::Orientation GetForient(const int i) const
Returns the orientation of face i with respect to the ordering of faces in the standard element.
Definition Geometry.h:395
Nektar::SpatialDomains::PointGeom::Sub
void Sub(PointGeom &a, PointGeom &b)
Definition PointGeom.cpp:118
Nektar::SpatialDomains::PointGeom::Mult
void Mult(PointGeom &a, PointGeom &b)
_this = a x b
Definition PointGeom.cpp:127
Nektar::SpatialDomains::PointGeom::dot
NekDouble dot(PointGeom &a)
retun the dot product between this and input a
Definition PointGeom.cpp:181
Nektar::SpatialDomains::PointGeom::UpdatePosition
void UpdatePosition(NekDouble x, NekDouble y, NekDouble z)
Definition PointGeom.cpp:101
Nektar::StdRegions::StdExpansion3D::GetEdgeNcoeffs
int GetEdgeNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th edge.
Definition StdExpansion3D.h:136
Nektar::StdRegions::StdExpansion3D::GetNedges
int GetNedges() const
return the number of edges in 3D expansion
Definition StdExpansion3D.h:121
Nektar::StdRegions::StdExpansion3D::GetEdgeInteriorToElementMap
void GetEdgeInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation traceOrient=eForwards)
Definition StdExpansion3D.h:141
Nektar::StdRegions::StdExpansion::GetBoundaryMap
void GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition StdExpansion.h:678
Nektar::StdRegions::StdExpansion::GetNcoeffs
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition StdExpansion.h:124
Nektar::StdRegions::StdExpansion::GetTotPoints
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition StdExpansion.h:134
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:500
Nektar::StdRegions::StdExpansion::DropLocStaticCondMatrix
void DropLocStaticCondMatrix(const LocalRegions::MatrixKey &mkey)
Definition StdExpansion.h:656
Nektar::StdRegions::StdExpansion::GetTraceNumPoints
int GetTraceNumPoints(const int i) const
This function returns the number of quadrature points belonging to the i-th trace.
Definition StdExpansion.h:281
Nektar::StdRegions::StdExpansion::GetStdMatrix
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition StdExpansion.h:612
Nektar::StdRegions::StdExpansion::GetPointsKeys
const LibUtilities::PointsKeyVector GetPointsKeys() const
Definition StdExpansion.h:1108
Nektar::StdRegions::StdExpansion::GetTraceIntNcoeffs
int GetTraceIntNcoeffs(const int i) const
Definition StdExpansion.h:266
Nektar::StdRegions::StdExpansion::GetLocStaticCondMatrix
DNekScalBlkMatSharedPtr GetLocStaticCondMatrix(const LocalRegions::MatrixKey &mkey)
Definition StdExpansion.h:650
Nektar::StdRegions::StdExpansion::GetVertexMap
int GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
Definition StdExpansion.h:688
Nektar::StdRegions::StdExpansion::NormVectorIProductWRTBase
void NormVectorIProductWRTBase(const Array< OneD, const NekDouble > &Fx, Array< OneD, NekDouble > &outarray)
Definition StdExpansion.h:622
Nektar::StdRegions::StdExpansion::GetLinStdExp
std::shared_ptr< StdExpansion > GetLinStdExp(void) const
Definition StdExpansion.h:380
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:537
Nektar::StdRegions::StdExpansion::GetNtraces
int GetNtraces() const
Returns the number of trace elements connected to this element.
Definition StdExpansion.h:354
Nektar::StdRegions::StdExpansion::GetNverts
int GetNverts() const
This function returns the number of vertices of the expansion domain.
Definition StdExpansion.h:246
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:693
Nektar::StdRegions::StdExpansion::DetShapeType
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition StdExpansion.h:370
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:717
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:261
Nektar::StdRegions::StdExpansion::GenMatrix
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
Definition StdExpansion.h:853
Nektar::StdRegions::StdExpansion::GetTraceNumModes
void GetTraceNumModes(const int tid, int &numModes0, int &numModes1, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2)
Definition StdExpansion.h:725
Nektar::StdRegions::StdExpansion::FwdTrans
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space.
Definition StdExpansion.h:1814
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:301
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:169
Nektar::StdRegions::StdExpansion::PhysDeriv
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
Definition StdExpansion.h:858
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition StdExpansion.h:1186
Nektar::StdRegions::StdMatrixKey::GetShapeType
LibUtilities::ShapeType GetShapeType() const
Definition StdMatrixKey.h:88
Nektar::StdRegions::StdMatrixKey::GetVarCoeffs
const VarCoeffMap & GetVarCoeffs() const
Definition StdMatrixKey.h:169
Nektar::StdRegions::StdMatrixKey::GetMatrixType
MatrixType GetMatrixType() const
Definition StdMatrixKey.h:83
Nektar::StdRegions::StdMatrixKey::HasVarCoeff
bool HasVarCoeff(const StdRegions::VarCoeffType &coeff) const
Definition StdMatrixKey.h:174
Nektar::StdRegions::StdMatrixKey::GetConstFactors
const ConstFactorMap & GetConstFactors() const
Definition StdMatrixKey.h:138
Nektar::StdRegions::StdMatrixKey::GetVarCoeff
const Array< OneD, const NekDouble > & GetVarCoeff(const StdRegions::VarCoeffType &coeff) const
Definition StdMatrixKey.h:148
Nektar::StdRegions::StdMatrixKey::GetConstFactor
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition StdMatrixKey.h:124
Nektar::StdRegions::StdMatrixKey::ConstFactorExists
bool ConstFactorExists(const ConstFactorType &factor) const
Definition StdMatrixKey.h:133
Nektar::LibUtilities::StdQuadData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb)
Definition ShapeType.hpp:133
Nektar::LibUtilities::StdTriData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb)
Definition ShapeType.hpp:109
Nektar::LibUtilities::Interp2D
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis,...
Definition Interp.cpp:101
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition Points.h:231
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::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::StdExpansionSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition StdExpansion.h:1783
Nektar::StdRegions::NullConstFactorMap
static ConstFactorMap NullConstFactorMap
Definition StdRegions.hpp:431
Nektar::StdRegions::NullVarCoeffMap
static VarCoeffMap NullVarCoeffMap
Definition StdRegions.hpp:376
Nektar::StdRegions::VarCoeffMap
std::map< StdRegions::VarCoeffType, VarCoeffEntry > VarCoeffMap
Definition StdRegions.hpp:375
Nektar::DNekScalMatSharedPtr
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Definition NekMatrixFwd.hpp:66
Nektar::DNekScalBlkMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition NekTypeDefs.hpp:79
Nektar::NullDNekMatSharedPtr
static DNekMatSharedPtr NullDNekMatSharedPtr
Definition NekTypeDefs.hpp:83
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::Gathr
void Gathr(I n, const T *x, const I *y, T *z)
Gather vector z[i] = x[y[i]].
Definition Vmath.hpp:507
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::Scatr
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition Vmath.hpp:539
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::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition Vmath.hpp:825
std
STL namespace.
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition scalar.hpp:290
Generated on Sun Nov 9 2025 20:16:58 for Nektar++ by doxygen 1.9.8