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