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