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