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