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