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