Nektar++
Expansion2D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File Expansion2D.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 Expansion2D routines
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <boost/core/ignore_unused.hpp>
36 
43 #include <LocalRegions/MatrixKey.h>
44 #include <LocalRegions/SegExp.h>
47 
48 using namespace std;
49 
50 namespace Nektar
51 {
52 namespace LocalRegions
53 {
54 Expansion2D::Expansion2D(SpatialDomains::Geometry2DSharedPtr pGeom)
55  : StdExpansion(), Expansion(pGeom), StdExpansion2D()
56 {
57 }
58 
60 {
61  DNekScalMatSharedPtr returnval;
63 
65  "Geometric information is not set up");
66 
67  switch (mkey.GetMatrixType())
68  {
69  case StdRegions::eMass:
70  {
71  if ((m_metricinfo->GetGtype() == SpatialDomains::eDeformed) ||
72  (mkey.GetNVarCoeff()))
73  {
74  NekDouble one = 1.0;
75  DNekMatSharedPtr mat = GenMatrix(mkey);
76 
77  returnval =
79  }
80  else
81  {
82  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
83  DNekMatSharedPtr mat = GetStdMatrix(mkey);
84 
85  returnval =
87  }
88  }
89  break;
91  {
92  MatrixKey masskey(mkey, StdRegions::eMass);
93  DNekScalMat &MassMat = *GetLocMatrix(masskey);
94 
95  // Generate a local copy of traceMat
98 
100  "Need to specify eFactorGJP to construct "
101  "a HelmholtzGJP matrix");
102 
104 
105  factor /= MassMat.Scale();
106 
107  int ntot = MassMat.GetRows() * MassMat.GetColumns();
108 
109  Vmath::Svtvp(ntot, factor, &NDTraceMat->GetPtr()[0], 1,
110  MassMat.GetRawPtr(), 1, &NDTraceMat->GetPtr()[0], 1);
111 
113  MassMat.Scale(), NDTraceMat);
114  }
115  break;
117  {
118  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
119  {
120  NekDouble one = 1.0;
122  DetShapeType(), *this);
123  DNekMatSharedPtr mat = GenMatrix(masskey);
124  mat->Invert();
125 
126  returnval =
128  }
129  else
130  {
131  NekDouble fac = 1.0 / (m_metricinfo->GetJac(ptsKeys))[0];
132  DNekMatSharedPtr mat = GetStdMatrix(mkey);
133 
134  returnval =
136  }
137  }
138  break;
142  {
143  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
144  mkey.GetNVarCoeff())
145  {
146  NekDouble one = 1.0;
147  DNekMatSharedPtr mat = GenMatrix(mkey);
148 
149  returnval =
151  }
152  else
153  {
154  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
156  m_metricinfo->GetDerivFactors(ptsKeys);
157  int dir = 0;
159  dir = 0;
161  dir = 1;
163  dir = 2;
164 
166  mkey.GetShapeType(), *this);
168  mkey.GetShapeType(), *this);
169 
170  DNekMat &deriv0 = *GetStdMatrix(deriv0key);
171  DNekMat &deriv1 = *GetStdMatrix(deriv1key);
172 
173  int rows = deriv0.GetRows();
174  int cols = deriv1.GetColumns();
175 
176  DNekMatSharedPtr WeakDeriv =
178  (*WeakDeriv) =
179  df[2 * dir][0] * deriv0 + df[2 * dir + 1][0] * deriv1;
180 
182  jac, WeakDeriv);
183  }
184  }
185  break;
187  {
188  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
189  mkey.GetNVarCoeff())
190  {
191  NekDouble one = 1.0;
192  DNekMatSharedPtr mat = GenMatrix(mkey);
193 
194  returnval =
196  }
197  else
198  {
199  int shapedim = 2;
200 
201  // dfdirxi = tan_{xi_x} * d \xi/dx
202  // + tan_{xi_y} * d \xi/dy
203  // + tan_{xi_z} * d \xi/dz
204  // dfdireta = tan_{eta_x} * d \eta/dx
205  // + tan_{xi_y} * d \xi/dy
206  // + tan_{xi_z} * d \xi/dz
207  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
209  m_metricinfo->GetDerivFactors(ptsKeys);
210 
211  Array<OneD, NekDouble> direction =
213 
214  // d / dx = df[0]*deriv0 + df[1]*deriv1
215  // d / dy = df[2]*deriv0 + df[3]*deriv1
216  // d / dz = df[4]*deriv0 + df[5]*deriv1
217 
218  // dfdir[dir] = e \cdot (d/dx, d/dy, d/dz)
219  // = (e^0 * df[0] + e^1 * df[2]
220  // + e^2 * df[4]) * deriv0
221  // + (e^0 * df[1] + e^1 * df[3]
222  // + e^2 * df[5]) * deriv1
223  // dfdir[dir] = e^0 * df[2 * 0 + dir]
224  // + e^1 * df[2 * 1 + dir]
225  // + e^2 * df [ 2 * 2 + dir]
226  Array<OneD, Array<OneD, NekDouble>> dfdir(shapedim);
227  Expansion::ComputeGmatcdotMF(df, direction, dfdir);
228 
229  StdRegions::VarCoeffMap dfdirxi;
230  StdRegions::VarCoeffMap dfdireta;
231 
232  dfdirxi[StdRegions::eVarCoeffWeakDeriv] = dfdir[0];
233  dfdireta[StdRegions::eVarCoeffWeakDeriv] = dfdir[1];
234 
236  mkey.GetShapeType(), *this,
239  mkey.GetShapeType(), *this,
241 
242  DNekMat &derivxi = *GetStdMatrix(derivxikey);
243  DNekMat &deriveta = *GetStdMatrix(derivetakey);
244 
245  int rows = derivxi.GetRows();
246  int cols = deriveta.GetColumns();
247 
248  DNekMatSharedPtr WeakDirDeriv =
250 
251  (*WeakDirDeriv) = derivxi + deriveta;
252 
253  // Add (\nabla \cdot e^k ) Mass
254  StdRegions::VarCoeffMap DiveMass;
255  DiveMass[StdRegions::eVarCoeffMass] =
257  StdRegions::StdMatrixKey stdmasskey(
258  StdRegions::eMass, mkey.GetShapeType(), *this,
260 
261  DNekMatSharedPtr DiveMassmat = GetStdMatrix(stdmasskey);
262 
263  (*WeakDirDeriv) = (*WeakDirDeriv) + (*DiveMassmat);
264 
266  jac, WeakDirDeriv);
267  }
268  break;
269  }
271  {
272  if ((m_metricinfo->GetGtype() == SpatialDomains::eDeformed) ||
273  (mkey.GetNVarCoeff() > 0) ||
275  {
276  NekDouble one = 1.0;
277  DNekMatSharedPtr mat = GenMatrix(mkey);
278 
279  returnval =
281  }
282  else
283  {
285  mkey.GetShapeType(), *this);
287  mkey.GetShapeType(), *this);
289  mkey.GetShapeType(), *this);
290 
291  DNekMat &lap00 = *GetStdMatrix(lap00key);
292  DNekMat &lap01 = *GetStdMatrix(lap01key);
293  DNekMat &lap11 = *GetStdMatrix(lap11key);
294 
295  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
297  m_metricinfo->GetGmat(ptsKeys);
298 
299  int rows = lap00.GetRows();
300  int cols = lap00.GetColumns();
301 
302  DNekMatSharedPtr lap =
304 
305  (*lap) = gmat[0][0] * lap00 +
306  gmat[1][0] * (lap01 + Transpose(lap01)) +
307  gmat[3][0] * lap11;
308 
309  returnval =
311  }
312  }
313  break;
315  {
316  DNekMatSharedPtr mat = GenMatrix(mkey);
317 
318  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, mat);
319  }
320  break;
322  {
324 
325  MatrixKey masskey(mkey, StdRegions::eMass);
326  DNekScalMat &MassMat = *GetLocMatrix(masskey);
327 
328  MatrixKey lapkey(mkey, StdRegions::eLaplacian);
329  DNekScalMat &LapMat = *GetLocMatrix(lapkey);
330 
331  int rows = LapMat.GetRows();
332  int cols = LapMat.GetColumns();
333 
334  DNekMatSharedPtr helm =
336 
337  NekDouble one = 1.0;
338  (*helm) = LapMat + factor * MassMat;
339 
340  returnval =
342  }
343  break;
345  {
346  MatrixKey helmkey(mkey, StdRegions::eHelmholtz);
347  DNekScalMat &HelmMat = *GetLocMatrix(helmkey);
348 
349  // Generate a local copy of traceMat
351  DNekMatSharedPtr NDTraceMat = Expansion2D::v_GenMatrix(key);
352 
354  "Need to specify eFactorGJP to construct "
355  "a HelmholtzGJP matrix");
356 
358 
359  factor /= HelmMat.Scale();
360 
361  int ntot = HelmMat.GetRows() * HelmMat.GetColumns();
362 
363  Vmath::Svtvp(ntot, factor, &NDTraceMat->GetPtr()[0], 1,
364  HelmMat.GetRawPtr(), 1, &NDTraceMat->GetPtr()[0], 1);
365 
367  HelmMat.Scale(), NDTraceMat);
368  }
369  break;
371  {
372  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
373  {
374  NekDouble one = 1.0;
375  DNekMatSharedPtr mat = GenMatrix(mkey);
376 
377  returnval =
379  }
380  else
381  {
382  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
383  DNekMatSharedPtr mat = GetStdMatrix(mkey);
384 
385  returnval =
387  }
388  }
389  break;
393  {
394  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
395  {
396  NekDouble one = 1.0;
397  DNekMatSharedPtr mat = GenMatrix(mkey);
398 
399  returnval =
401  }
402  else
403  {
404  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
405 
406  const Array<TwoD, const NekDouble> &df =
407  m_metricinfo->GetDerivFactors(ptsKeys);
408  int dir = 0;
410  dir = 0;
412  dir = 1;
414  dir = 2;
415 
417  mkey.GetShapeType(), *this);
419  mkey.GetShapeType(), *this);
420 
421  DNekMat &stdiprod0 = *GetStdMatrix(iProdDeriv0Key);
422  DNekMat &stdiprod1 = *GetStdMatrix(iProdDeriv0Key);
423 
424  int rows = stdiprod0.GetRows();
425  int cols = stdiprod1.GetColumns();
426 
427  DNekMatSharedPtr mat =
429  (*mat) =
430  df[2 * dir][0] * stdiprod0 + df[2 * dir + 1][0] * stdiprod1;
431 
432  returnval =
434  }
435  }
436  break;
437 
439  {
440  NekDouble one = 1.0;
441 
443  *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
444 
445  DNekMatSharedPtr mat = GenMatrix(hkey);
446 
447  mat->Invert();
448 
449  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
450  }
451  break;
453  {
455  "Matrix only setup for quad elements currently");
456  DNekMatSharedPtr m_Ix;
457  Array<OneD, NekDouble> coords(1, 0.0);
459  int edge = static_cast<int>(factors[StdRegions::eFactorGaussEdge]);
460 
461  coords[0] = (edge == 0 || edge == 3) ? -1.0 : 1.0;
462 
463  m_Ix = m_base[(edge + 1) % 2]->GetI(coords);
464 
465  returnval =
467  }
468  break;
470  {
471  NekDouble one = 1.0;
473  *this, mkey.GetConstFactors(),
474  mkey.GetVarCoeffs());
475  DNekScalBlkMatSharedPtr helmStatCond =
476  GetLocStaticCondMatrix(helmkey);
477  DNekScalMatSharedPtr A = helmStatCond->GetBlock(0, 0);
479 
481  }
482  break;
483  default:
484  {
485  NekDouble one = 1.0;
486  DNekMatSharedPtr mat = GenMatrix(mkey);
487 
488  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
489  }
490  break;
491  }
492 
493  return returnval;
494 }
495 
497  const int edge, const ExpansionSharedPtr &EdgeExp,
500 {
501  ASSERTL1(GetCoordim() == 2, "Routine only set up for two-dimensions");
502 
504  GetTraceNormal(edge);
505 
506  if (m_requireNeg.size() == 0)
507  {
508  int nedges = GetNtraces();
509  m_requireNeg.resize(nedges);
510 
511  for (int i = 0; i < nedges; ++i)
512  {
513  m_requireNeg[i] = false;
514 
515  ExpansionSharedPtr edgeExp = m_traceExp[i].lock();
516 
517  if (edgeExp->GetRightAdjacentElementExp())
518  {
519  if (edgeExp->GetRightAdjacentElementExp()
520  ->GetGeom()
521  ->GetGlobalID() == GetGeom()->GetGlobalID())
522  {
523  m_requireNeg[i] = true;
524  }
525  }
526  }
527  }
528 
529  // We allow the case of mixed polynomial order by supporting only
530  // those modes on the edge common to both adjoining elements. This
531  // is enforced here by taking the minimum size and padding with
532  // zeros.
533  int nquad_e = min(EdgeExp->GetNumPoints(0), int(normals[0].size()));
534 
535  int nEdgePts = EdgeExp->GetTotPoints();
536  Array<OneD, NekDouble> edgePhys(nEdgePts);
537  Vmath::Vmul(nquad_e, normals[0], 1, Fx, 1, edgePhys, 1);
538  Vmath::Vvtvp(nquad_e, normals[1], 1, Fy, 1, edgePhys, 1, edgePhys, 1);
539 
540  Expansion1DSharedPtr locExp = EdgeExp->as<Expansion1D>();
541 
542  if (m_requireNeg[edge])
543  {
544  if (locExp->GetRightAdjacentElementExp()->GetGeom()->GetGlobalID() ==
545  m_geom->GetGlobalID())
546  {
547  Vmath::Neg(nquad_e, edgePhys, 1);
548  }
549  }
550 
551  AddEdgeNormBoundaryInt(edge, EdgeExp, edgePhys, outarray);
552 }
553 
555  const int edge, const ExpansionSharedPtr &EdgeExp,
557 {
558  int i;
559 
560  if (m_requireNeg.size() == 0)
561  {
562  int nedges = GetNtraces();
563  m_requireNeg.resize(nedges);
564 
565  for (i = 0; i < nedges; ++i)
566  {
567  m_requireNeg[i] = false;
568 
569  ExpansionSharedPtr edgeExp = m_traceExp[i].lock();
570 
571  if (edgeExp->GetRightAdjacentElementExp())
572  {
573  if (edgeExp->GetRightAdjacentElementExp()
574  ->GetGeom()
575  ->GetGlobalID() == GetGeom()->GetGlobalID())
576  {
577  m_requireNeg[i] = true;
578  }
579  }
580  }
581  }
582 
584  GetBasisNumModes(1), 0, edge, GetTraceOrient(edge));
585 
587 
588  // Order of the element
589  int order_e = map->size();
590  // Order of the trace
591  int n_coeffs = EdgeExp->GetNcoeffs();
592 
593  Array<OneD, NekDouble> edgeCoeffs(n_coeffs);
594  if (n_coeffs != order_e) // Going to orthogonal space
595  {
596  EdgeExp->FwdTrans(Fn, edgeCoeffs);
597  Expansion1DSharedPtr locExp = EdgeExp->as<Expansion1D>();
598 
599  if (m_requireNeg[edge])
600  {
601  Vmath::Neg(n_coeffs, edgeCoeffs, 1);
602  }
603 
604  Array<OneD, NekDouble> coeff(n_coeffs, 0.0);
606  ((LibUtilities::BasisType)1); // 1-->Ortho_A
607  LibUtilities::BasisKey bkey_ortho(btype,
608  EdgeExp->GetBasis(0)->GetNumModes(),
609  EdgeExp->GetBasis(0)->GetPointsKey());
610  LibUtilities::BasisKey bkey(EdgeExp->GetBasis(0)->GetBasisType(),
611  EdgeExp->GetBasis(0)->GetNumModes(),
612  EdgeExp->GetBasis(0)->GetPointsKey());
613  LibUtilities::InterpCoeff1D(bkey, edgeCoeffs, bkey_ortho, coeff);
614 
615  // Cutting high frequencies
616  for (i = order_e; i < n_coeffs; i++)
617  {
618  coeff[i] = 0.0;
619  }
620 
621  LibUtilities::InterpCoeff1D(bkey_ortho, coeff, bkey, edgeCoeffs);
622 
624  LibUtilities::eSegment, *EdgeExp);
625  EdgeExp->MassMatrixOp(edgeCoeffs, edgeCoeffs, masskey);
626  }
627  else
628  {
629  EdgeExp->IProductWRTBase(Fn, edgeCoeffs);
630 
631  Expansion1DSharedPtr locExp = EdgeExp->as<Expansion1D>();
632 
633  if (m_requireNeg[edge])
634  {
635  Vmath::Neg(n_coeffs, edgeCoeffs, 1);
636  }
637  }
638 
639  // Implementation for all the basis except Gauss points
640  if (EdgeExp->GetBasis(0)->GetBasisType() != LibUtilities::eGauss_Lagrange)
641  {
642  // add data to outarray if forward edge normal is outwards
643  for (i = 0; i < order_e; ++i)
644  {
645  outarray[(*map)[i].index] += (*map)[i].sign * edgeCoeffs[i];
646  }
647  }
648  else
649  {
650  int nCoeffs0, nCoeffs1;
651  int j;
652 
654  factors[StdRegions::eFactorGaussEdge] = edge;
656  *this, factors);
657 
658  DNekMatSharedPtr mat_gauss = m_stdMatrixManager[key];
659 
660  switch (edge)
661  {
662  case 0:
663  {
664  nCoeffs1 = m_base[1]->GetNumModes();
665 
666  for (i = 0; i < order_e; ++i)
667  {
668  for (j = 0; j < nCoeffs1; j++)
669  {
670  outarray[(*map)[i].index + j * order_e] +=
671  mat_gauss->GetPtr()[j] * (*map)[i].sign *
672  edgeCoeffs[i];
673  }
674  }
675  break;
676  }
677  case 1:
678  {
679  nCoeffs0 = m_base[0]->GetNumModes();
680 
681  for (i = 0; i < order_e; ++i)
682  {
683  for (j = 0; j < nCoeffs0; j++)
684  {
685  outarray[(*map)[i].index - j] +=
686  mat_gauss->GetPtr()[order_e - 1 - j] *
687  (*map)[i].sign * edgeCoeffs[i];
688  }
689  }
690  break;
691  }
692  case 2:
693  {
694  nCoeffs1 = m_base[1]->GetNumModes();
695 
696  for (i = 0; i < order_e; ++i)
697  {
698  for (j = 0; j < nCoeffs1; j++)
699  {
700  outarray[(*map)[i].index - j * order_e] +=
701  mat_gauss->GetPtr()[order_e - 1 - j] *
702  (*map)[i].sign * edgeCoeffs[i];
703  }
704  }
705  break;
706  }
707  case 3:
708  {
709  nCoeffs0 = m_base[0]->GetNumModes();
710 
711  for (i = 0; i < order_e; ++i)
712  {
713  for (j = 0; j < nCoeffs0; j++)
714  {
715  outarray[(*map)[i].index + j] +=
716  mat_gauss->GetPtr()[j] * (*map)[i].sign *
717  edgeCoeffs[i];
718  }
719  }
720  break;
721  }
722  default:
723  ASSERTL0(false, "edge value (< 3) is out of range");
724  break;
725  }
726  }
727 }
728 
731 {
732  int i, cnt = 0;
733  int nedges = GetNtraces();
735 
736  for (i = 0; i < nedges; ++i)
737  {
738  EdgeExp[i]->SetCoeffsToOrientation(
739  GetTraceOrient(i), e_tmp = inout + cnt, e_tmp = inout + cnt);
740  cnt += GetTraceNcoeffs(i);
741  }
742 }
743 
744 /**
745  * Computes the C matrix entries due to the presence of the identity
746  * matrix in Eqn. 32.
747  */
748 void Expansion2D::AddNormTraceInt(const int dir,
751  Array<OneD, NekDouble> &outarray,
752  const StdRegions::VarCoeffMap &varcoeffs)
753 {
754  int i, e, cnt;
755  int order_e, nquad_e;
756  int nedges = GetNtraces();
757 
758  cnt = 0;
759  for (e = 0; e < nedges; ++e)
760  {
761  order_e = EdgeExp[e]->GetNcoeffs();
762  nquad_e = EdgeExp[e]->GetNumPoints(0);
763 
765  GetTraceNormal(e);
766  Array<OneD, NekDouble> edgeCoeffs(order_e);
767  Array<OneD, NekDouble> edgePhys(nquad_e);
768 
769  for (i = 0; i < order_e; ++i)
770  {
771  edgeCoeffs[i] = inarray[i + cnt];
772  }
773  cnt += order_e;
774 
775  EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
776 
777  // Multiply by variable coefficient
778  /// @TODO: Document this
779  // StdRegions::VarCoeffType VarCoeff[3] = {StdRegions::eVarCoeffD00,
780  // StdRegions::eVarCoeffD11,
781  // StdRegions::eVarCoeffD22};
782  // StdRegions::VarCoeffMap::const_iterator x;
783  // Array<OneD, NekDouble> varcoeff_work(nquad_e);
784 
785  // if ((x = varcoeffs.find(VarCoeff[dir])) != varcoeffs.end())
786  // {
787  // GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
788  // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
789  // }
790 
791  if (varcoeffs.find(StdRegions::eVarCoeffMF1x) != varcoeffs.end())
792  {
793  // MMF case
794  Array<OneD, NekDouble> ncdotMF_e =
795  v_GetnEdgecdotMF(dir, e, EdgeExp[e], normals, varcoeffs);
796 
797  Vmath::Vmul(nquad_e, ncdotMF_e, 1, edgePhys, 1, edgePhys, 1);
798  }
799  else
800  {
801  Vmath::Vmul(nquad_e, normals[dir], 1, edgePhys, 1, edgePhys, 1);
802  }
803 
804  AddEdgeBoundaryInt(e, EdgeExp[e], edgePhys, outarray, varcoeffs);
805  }
806 }
807 
809  const int dir, Array<OneD, ExpansionSharedPtr> &EdgeExp,
810  Array<OneD, Array<OneD, NekDouble>> &edgeCoeffs,
811  Array<OneD, NekDouble> &outarray)
812 {
813  int e;
814  int nquad_e;
815  int nedges = GetNtraces();
816 
817  for (e = 0; e < nedges; ++e)
818  {
819  nquad_e = EdgeExp[e]->GetNumPoints(0);
820 
821  Array<OneD, NekDouble> edgePhys(nquad_e);
823  GetTraceNormal(e);
824 
825  EdgeExp[e]->BwdTrans(edgeCoeffs[e], edgePhys);
826 
827  Vmath::Vmul(nquad_e, normals[dir], 1, edgePhys, 1, edgePhys, 1);
828 
829  AddEdgeBoundaryInt(e, EdgeExp[e], edgePhys, outarray);
830  }
831 }
832 
833 /**
834  * For a given edge add the \tilde{F}_1j contributions
835  */
837  ExpansionSharedPtr &EdgeExp,
838  Array<OneD, NekDouble> &edgePhys,
839  Array<OneD, NekDouble> &outarray,
840  const StdRegions::VarCoeffMap &varcoeffs)
841 {
842  int i;
843  int order_e = EdgeExp->GetNcoeffs();
844  int nquad_e = EdgeExp->GetNumPoints(0);
847  Array<OneD, NekDouble> coeff(order_e);
848 
849  GetTraceToElementMap(edge, map, sign, v_GetTraceOrient(edge));
850 
854  StdRegions::VarCoeffMap::const_iterator x;
855 
856  /// @TODO Variable coeffs
857  if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
858  {
859  Array<OneD, NekDouble> work(nquad_e);
860  GetPhysEdgeVarCoeffsFromElement(edge, EdgeExp, x->second, work);
861  Vmath::Vmul(nquad_e, work, 1, edgePhys, 1, edgePhys, 1);
862  }
863 
864  EdgeExp->IProductWRTBase(edgePhys, coeff);
865 
866  // add data to out array
867  for (i = 0; i < order_e; ++i)
868  {
869  outarray[map[i]] += sign[i] * coeff[i];
870  }
871 }
872 
873 // This method assumes that data in EdgeExp is orientated according to
874 // elemental counter clockwise format AddHDGHelmholtzTraceTerms with
875 // directions
877  const NekDouble tau, const Array<OneD, const NekDouble> &inarray,
879  const StdRegions::VarCoeffMap &dirForcing, Array<OneD, NekDouble> &outarray)
880 {
881  ASSERTL0(&inarray[0] != &outarray[0],
882  "Input and output arrays use the same memory");
883 
884  int e, cnt, order_e, nedges = GetNtraces();
886 
887  cnt = 0;
888 
889  for (e = 0; e < nedges; ++e)
890  {
891  order_e = EdgeExp[e]->GetNcoeffs();
892  Array<OneD, NekDouble> edgeCoeffs(order_e);
893  Array<OneD, NekDouble> edgePhys(EdgeExp[e]->GetTotPoints());
894 
895  Vmath::Vcopy(order_e, tmp = inarray + cnt, 1, edgeCoeffs, 1);
896  EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
897  AddHDGHelmholtzEdgeTerms(tau, e, EdgeExp, edgePhys, dirForcing,
898  outarray);
899 
900  cnt += order_e;
901  }
902 }
903 
904 // evaluate additional terms in HDG edges. Not that this assumes that
905 // edges are unpacked into local cartesian order.
907  const NekDouble tau, const int edge,
909  const StdRegions::VarCoeffMap &varcoeffs, Array<OneD, NekDouble> &outarray)
910 {
911  bool mmf = (varcoeffs.find(StdRegions::eVarCoeffMF1x) != varcoeffs.end());
912  int i, j, n;
913  int nquad_e = EdgeExp[edge]->GetNumPoints(0);
914  int order_e = EdgeExp[edge]->GetNcoeffs();
915  int coordim = mmf ? 2 : GetCoordim();
916  int ncoeffs = GetNcoeffs();
917 
918  Array<OneD, NekDouble> inval(nquad_e);
919  Array<OneD, NekDouble> outcoeff(order_e);
920  Array<OneD, NekDouble> tmpcoeff(ncoeffs);
921 
923  GetTraceNormal(edge);
924 
927 
928  StdRegions::Orientation edgedir = GetTraceOrient(edge);
929 
930  DNekVec Coeffs(ncoeffs, outarray, eWrapper);
931  DNekVec Tmpcoeff(ncoeffs, tmpcoeff, eWrapper);
932 
933  GetTraceToElementMap(edge, emap, sign, edgedir);
934 
938 
942 
943  StdRegions::VarCoeffMap::const_iterator x;
944  /// @TODO: What direction to use here??
945  if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
946  {
947  Array<OneD, NekDouble> work(nquad_e);
948  GetPhysEdgeVarCoeffsFromElement(edge, EdgeExp[edge], x->second, work);
949  Vmath::Vmul(nquad_e, work, 1, edgePhys, 1, edgePhys, 1);
950  }
951 
952  //================================================================
953  // Add F = \tau <phi_i,in_phys>
954  // Fill edge and take inner product
955  EdgeExp[edge]->IProductWRTBase(edgePhys, outcoeff);
956  // add data to out array
957  for (i = 0; i < order_e; ++i)
958  {
959  outarray[emap[i]] += sign[i] * tau * outcoeff[i];
960  }
961  //================================================================
962 
963  //===============================================================
964  // Add -\sum_i D_i^T M^{-1} G_i + E_i M^{-1} G_i =
965  // \sum_i D_i M^{-1} G_i term
966 
967  // Two independent direction
968  DNekScalMatSharedPtr invMass;
969  for (n = 0; n < coordim; ++n)
970  {
971  if (mmf)
972  {
974  Weight[StdRegions::eVarCoeffMass] = v_GetMFMag(n, varcoeffs);
975 
976  MatrixKey invMasskey(StdRegions::eInvMass, DetShapeType(), *this,
978 
979  invMass = GetLocMatrix(invMasskey);
980 
981  Array<OneD, NekDouble> ncdotMF_e =
982  v_GetnEdgecdotMF(n, edge, EdgeExp[edge], normals, varcoeffs);
983 
984  Vmath::Vmul(nquad_e, ncdotMF_e, 1, edgePhys, 1, inval, 1);
985  }
986  else
987  {
988  Vmath::Vmul(nquad_e, normals[n], 1, edgePhys, 1, inval, 1);
990  }
991 
992  // Multiply by variable coefficient
993  /// @TODO: Document this (probably not needed)
994  // StdRegions::VarCoeffMap::const_iterator x;
995  // if ((x = varcoeffs.find(VarCoeff[n])) !=
996  // varcoeffs.end())
997  // {
998  // GetPhysEdgeVarCoeffsFromElement(edge,EdgeExp[edge],x->second,varcoeff_work);
999  // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[edge]->GetPhys(),1,EdgeExp[edge]->UpdatePhys(),1);
1000  // }
1001 
1002  EdgeExp[edge]->IProductWRTBase(inval, outcoeff);
1003 
1004  // M^{-1} G
1005  for (i = 0; i < ncoeffs; ++i)
1006  {
1007  tmpcoeff[i] = 0;
1008  for (j = 0; j < order_e; ++j)
1009  {
1010  tmpcoeff[i] += (*invMass)(i, emap[j]) * sign[j] * outcoeff[j];
1011  }
1012  }
1013 
1014  if (mmf)
1015  {
1016  StdRegions::VarCoeffMap VarCoeffDirDeriv;
1017  VarCoeffDirDeriv[StdRegions::eVarCoeffMF] =
1018  v_GetMF(n, coordim, varcoeffs);
1019  VarCoeffDirDeriv[StdRegions::eVarCoeffMFDiv] =
1020  v_GetMFDiv(n, varcoeffs);
1021 
1024  VarCoeffDirDeriv);
1025 
1026  DNekScalMat &Dmat = *GetLocMatrix(Dmatkey);
1027 
1028  Coeffs = Coeffs + Dmat * Tmpcoeff;
1029  }
1030  else
1031  {
1032  if (varcoeffs.find(VarCoeff[n]) != varcoeffs.end())
1033  {
1034  MatrixKey mkey(DerivType[n], DetShapeType(), *this,
1035  StdRegions::NullConstFactorMap, varcoeffs);
1036 
1037  DNekScalMat &Dmat = *GetLocMatrix(mkey);
1038  Coeffs = Coeffs + Dmat * Tmpcoeff;
1039  }
1040  else
1041  {
1042  DNekScalMat &Dmat = *GetLocMatrix(DerivType[n]);
1043  Coeffs = Coeffs + Dmat * Tmpcoeff;
1044  }
1045  }
1046  }
1047 }
1048 
1049 /**
1050  * Extracts the variable coefficients along an edge
1051  */
1053  const int edge, ExpansionSharedPtr &EdgeExp,
1054  const Array<OneD, const NekDouble> &varcoeff,
1055  Array<OneD, NekDouble> &outarray)
1056 {
1058  Array<OneD, NekDouble> edgetmp(EdgeExp->GetNcoeffs());
1059 
1060  // FwdTrans varcoeffs
1061  FwdTrans(varcoeff, tmp);
1062 
1063  // Map to edge
1066  StdRegions::Orientation edgedir = GetTraceOrient(edge);
1067  GetTraceToElementMap(edge, emap, sign, edgedir);
1068 
1069  for (unsigned int i = 0; i < EdgeExp->GetNcoeffs(); ++i)
1070  {
1071  edgetmp[i] = tmp[emap[i]];
1072  }
1073 
1074  // BwdTrans
1075  EdgeExp->BwdTrans(edgetmp, outarray);
1076 }
1077 
1078 /**
1079  * Computes matrices needed for the HDG formulation. References to
1080  * equations relate to the following paper:
1081  * R. M. Kirby, S. J. Sherwin, B. Cockburn, To CG or to HDG: A
1082  * Comparative Study, J. Sci. Comp P1-30
1083  * DOI 10.1007/s10915-011-9501-7
1084  */
1086 {
1087  DNekMatSharedPtr returnval;
1088 
1089  switch (mkey.GetMatrixType())
1090  {
1091  // (Z^e)^{-1} (Eqn. 33, P22)
1093  {
1095  "HybridDGHelmholtz matrix not set up "
1096  "for non boundary-interior expansions");
1097 
1098  int i, j, k;
1099  NekDouble lambdaval =
1102  int ncoeffs = GetNcoeffs();
1103  int nedges = GetNtraces();
1104  int shapedim = 2;
1105  const StdRegions::VarCoeffMap &varcoeffs = mkey.GetVarCoeffs();
1106  bool mmf =
1107  (varcoeffs.find(StdRegions::eVarCoeffMF1x) != varcoeffs.end());
1108 
1112  ExpansionSharedPtr EdgeExp;
1113 
1114  int order_e, coordim = GetCoordim();
1119  DNekMat LocMat(ncoeffs, ncoeffs);
1120 
1121  returnval =
1123  DNekMat &Mat = *returnval;
1124  Vmath::Zero(ncoeffs * ncoeffs, Mat.GetPtr(), 1);
1125 
1129 
1130  StdRegions::VarCoeffMap::const_iterator x;
1131 
1132  for (i = 0; i < coordim; ++i)
1133  {
1134  if (mmf)
1135  {
1136  if (i < shapedim)
1137  {
1138  StdRegions::VarCoeffMap VarCoeffDirDeriv;
1139  VarCoeffDirDeriv[StdRegions::eVarCoeffMF] =
1140  v_GetMF(i, shapedim, varcoeffs);
1141  VarCoeffDirDeriv[StdRegions::eVarCoeffMFDiv] =
1142  v_GetMFDiv(i, varcoeffs);
1143 
1145  DetShapeType(), *this,
1147  VarCoeffDirDeriv);
1148 
1149  DNekScalMat &Dmat = *GetLocMatrix(Dmatkey);
1150 
1151  StdRegions::VarCoeffMap Weight;
1152  Weight[StdRegions::eVarCoeffMass] =
1153  v_GetMFMag(i, mkey.GetVarCoeffs());
1154 
1155  MatrixKey invMasskey(
1158 
1159  DNekScalMat &invMass = *GetLocMatrix(invMasskey);
1160 
1161  Mat = Mat + Dmat * invMass * Transpose(Dmat);
1162  }
1163  }
1164  else if (mkey.HasVarCoeff(Coeffs[i]))
1165  {
1166  MatrixKey DmatkeyL(DerivType[i], DetShapeType(), *this,
1168  mkey.GetVarCoeffAsMap(Coeffs[i]));
1169 
1170  MatrixKey DmatkeyR(DerivType[i], DetShapeType(), *this);
1171 
1172  DNekScalMat &DmatL = *GetLocMatrix(DmatkeyL);
1173  DNekScalMat &DmatR = *GetLocMatrix(DmatkeyR);
1174  Mat = Mat + DmatL * invMass * Transpose(DmatR);
1175  }
1176  else
1177  {
1178  DNekScalMat &Dmat = *GetLocMatrix(DerivType[i]);
1179  Mat = Mat + Dmat * invMass * Transpose(Dmat);
1180  }
1181  }
1182 
1183  // Add Mass Matrix Contribution for Helmholtz problem
1185  Mat = Mat + lambdaval * Mass;
1186 
1187  // Add tau*E_l using elemental mass matrices on each edge
1188  for (i = 0; i < nedges; ++i)
1189  {
1190  EdgeExp = GetTraceExp(i);
1191  order_e = EdgeExp->GetNcoeffs();
1192 
1193  int nq = EdgeExp->GetNumPoints(0);
1194  GetTraceToElementMap(i, emap, sign, edgedir);
1195 
1196  // @TODO: Document
1197  StdRegions::VarCoeffMap edgeVarCoeffs;
1199  {
1200  Array<OneD, NekDouble> mu(nq);
1202  i, EdgeExp, mkey.GetVarCoeff(StdRegions::eVarCoeffD00),
1203  mu);
1204  edgeVarCoeffs[StdRegions::eVarCoeffMass] = mu;
1205  }
1206  DNekScalMat &eMass = *EdgeExp->GetLocMatrix(
1208  edgeVarCoeffs);
1209  // DNekScalMat &eMass =
1210  // *EdgeExp->GetLocMatrix(StdRegions::eMass);
1211 
1212  for (j = 0; j < order_e; ++j)
1213  {
1214  for (k = 0; k < order_e; ++k)
1215  {
1216  Mat(emap[j], emap[k]) =
1217  Mat(emap[j], emap[k]) +
1218  tau * sign[j] * sign[k] * eMass(j, k);
1219  }
1220  }
1221  }
1222  }
1223  break;
1224  // U^e (P22)
1226  {
1227  int i, j, k;
1228  int nbndry = NumDGBndryCoeffs();
1229  int ncoeffs = GetNcoeffs();
1230  int nedges = GetNtraces();
1232 
1233  Array<OneD, NekDouble> lambda(nbndry);
1234  DNekVec Lambda(nbndry, lambda, eWrapper);
1235  Array<OneD, NekDouble> ulam(ncoeffs);
1236  DNekVec Ulam(ncoeffs, ulam, eWrapper);
1237  Array<OneD, NekDouble> f(ncoeffs);
1238  DNekVec F(ncoeffs, f, eWrapper);
1239 
1240  Array<OneD, ExpansionSharedPtr> EdgeExp(nedges);
1241  // declare matrix space
1242  returnval =
1244  DNekMat &Umat = *returnval;
1245 
1246  // Z^e matrix
1248  *this, mkey.GetConstFactors(),
1249  mkey.GetVarCoeffs());
1250  DNekScalMat &invHmat = *GetLocMatrix(newkey);
1251 
1254 
1255  for (i = 0; i < nedges; ++i)
1256  {
1257  EdgeExp[i] = GetTraceExp(i);
1258  }
1259 
1260  // for each degree of freedom of the lambda space
1261  // calculate Umat entry
1262  // Generate Lambda to U_lambda matrix
1263  for (j = 0; j < nbndry; ++j)
1264  {
1265  // standard basis vectors e_j
1266  Vmath::Zero(nbndry, &lambda[0], 1);
1267  Vmath::Zero(ncoeffs, &f[0], 1);
1268  lambda[j] = 1.0;
1269 
1270  SetTraceToGeomOrientation(EdgeExp, lambda);
1271 
1272  // Compute F = [I D_1 M^{-1} D_2 M^{-1}] C e_j
1273  AddHDGHelmholtzTraceTerms(tau, lambda, EdgeExp,
1274  mkey.GetVarCoeffs(), f);
1275 
1276  // Compute U^e_j
1277  Ulam = invHmat * F; // generate Ulam from lambda
1278 
1279  // fill column of matrix
1280  for (k = 0; k < ncoeffs; ++k)
1281  {
1282  Umat(k, j) = Ulam[k];
1283  }
1284  }
1285  }
1286  break;
1287  // Q_0, Q_1, Q_2 matrices (P23)
1288  // Each are a product of a row of Eqn 32 with the C matrix.
1289  // Rather than explicitly computing all of Eqn 32, we note each
1290  // row is almost a multiple of U^e, so use that as our starting
1291  // point.
1295  {
1296  int i = 0;
1297  int j = 0;
1298  int k = 0;
1299  int dir = 0;
1300  int nbndry = NumDGBndryCoeffs();
1301  int ncoeffs = GetNcoeffs();
1302  int nedges = GetNtraces();
1303  int shapedim = 2;
1304 
1305  Array<OneD, NekDouble> lambda(nbndry);
1306  DNekVec Lambda(nbndry, lambda, eWrapper);
1307  Array<OneD, ExpansionSharedPtr> EdgeExp(nedges);
1308 
1309  Array<OneD, NekDouble> ulam(ncoeffs);
1310  DNekVec Ulam(ncoeffs, ulam, eWrapper);
1311  Array<OneD, NekDouble> f(ncoeffs);
1312  DNekVec F(ncoeffs, f, eWrapper);
1313 
1314  // declare matrix space
1315  returnval =
1317  DNekMat &Qmat = *returnval;
1318 
1319  // Lambda to U matrix
1321  *this, mkey.GetConstFactors(),
1322  mkey.GetVarCoeffs());
1323  DNekScalMat &lamToU = *GetLocMatrix(lamToUkey);
1324 
1325  // Inverse mass matrix
1327 
1328  for (i = 0; i < nedges; ++i)
1329  {
1330  EdgeExp[i] = GetTraceExp(i);
1331  }
1332 
1333  // Weak Derivative matrix
1334  DNekScalMatSharedPtr Dmat;
1335  switch (mkey.GetMatrixType())
1336  {
1338  dir = 0;
1340  break;
1342  dir = 1;
1344  break;
1346  dir = 2;
1348  break;
1349  default:
1350  ASSERTL0(false, "Direction not known");
1351  break;
1352  }
1353 
1354  const StdRegions::VarCoeffMap &varcoeffs = mkey.GetVarCoeffs();
1355  if (varcoeffs.find(StdRegions::eVarCoeffMF1x) != varcoeffs.end())
1356  {
1357  StdRegions::VarCoeffMap VarCoeffDirDeriv;
1358  VarCoeffDirDeriv[StdRegions::eVarCoeffMF] =
1359  v_GetMF(dir, shapedim, varcoeffs);
1360  VarCoeffDirDeriv[StdRegions::eVarCoeffMFDiv] =
1361  v_GetMFDiv(dir, varcoeffs);
1362 
1363  MatrixKey Dmatkey(
1365  StdRegions::NullConstFactorMap, VarCoeffDirDeriv);
1366 
1367  Dmat = GetLocMatrix(Dmatkey);
1368 
1369  StdRegions::VarCoeffMap Weight;
1370  Weight[StdRegions::eVarCoeffMass] =
1371  v_GetMFMag(dir, mkey.GetVarCoeffs());
1372 
1375  Weight);
1376 
1377  invMass = *GetLocMatrix(invMasskey);
1378  }
1379  else
1380  {
1384 
1385  Dmat = GetLocMatrix(DerivType[dir]);
1386 
1388  *this);
1389  invMass = *GetLocMatrix(invMasskey);
1390  }
1391 
1392  // for each degree of freedom of the lambda space
1393  // calculate Qmat entry
1394  // Generate Lambda to Q_lambda matrix
1395  for (j = 0; j < nbndry; ++j)
1396  {
1397  Vmath::Zero(nbndry, &lambda[0], 1);
1398  lambda[j] = 1.0;
1399 
1400  // for lambda[j] = 1 this is the solution to ulam
1401  for (k = 0; k < ncoeffs; ++k)
1402  {
1403  Ulam[k] = lamToU(k, j);
1404  }
1405 
1406  // -D^T ulam
1407  Vmath::Neg(ncoeffs, &ulam[0], 1);
1408  F = Transpose(*Dmat) * Ulam;
1409 
1410  SetTraceToGeomOrientation(EdgeExp, lambda);
1411 
1412  // Add the C terms resulting from the I's on the
1413  // diagonals of Eqn 32
1414  AddNormTraceInt(dir, lambda, EdgeExp, f, mkey.GetVarCoeffs());
1415 
1416  // finally multiply by inverse mass matrix
1417  Ulam = invMass * F;
1418 
1419  // fill column of matrix (Qmat is in column major format)
1420  Vmath::Vcopy(ncoeffs, &ulam[0], 1,
1421  &(Qmat.GetPtr())[0] + j * ncoeffs, 1);
1422  }
1423  }
1424  break;
1425  // Matrix K (P23)
1427  {
1428  int i, j, e, cnt;
1429  int order_e, nquad_e;
1430  int nbndry = NumDGBndryCoeffs();
1431  int coordim = GetCoordim();
1432  int nedges = GetNtraces();
1434  StdRegions::VarCoeffMap::const_iterator x;
1435  const StdRegions::VarCoeffMap &varcoeffs = mkey.GetVarCoeffs();
1436  bool mmf =
1437  (varcoeffs.find(StdRegions::eVarCoeffMF1x) != varcoeffs.end());
1438 
1439  Array<OneD, NekDouble> work, varcoeff_work;
1441  Array<OneD, ExpansionSharedPtr> EdgeExp(nedges);
1442  Array<OneD, NekDouble> lam(nbndry);
1443 
1446  StdRegions::Orientation edgedir;
1447 
1448  // declare matrix space
1449  returnval =
1451  DNekMat &BndMat = *returnval;
1452 
1453  DNekScalMatSharedPtr LamToQ[3];
1454 
1455  // Matrix to map Lambda to U
1457  *this, mkey.GetConstFactors(),
1458  mkey.GetVarCoeffs());
1459  DNekScalMat &LamToU = *GetLocMatrix(LamToUkey);
1460 
1461  // Matrix to map Lambda to Q0
1463  *this, mkey.GetConstFactors(),
1464  mkey.GetVarCoeffs());
1465  LamToQ[0] = GetLocMatrix(LamToQ0key);
1466 
1467  // Matrix to map Lambda to Q1
1469  *this, mkey.GetConstFactors(),
1470  mkey.GetVarCoeffs());
1471  LamToQ[1] = GetLocMatrix(LamToQ1key);
1472 
1473  // Matrix to map Lambda to Q2 for 3D coordinates
1474  if (coordim == 3)
1475  {
1476  MatrixKey LamToQ2key(
1478  mkey.GetConstFactors(), mkey.GetVarCoeffs());
1479  LamToQ[2] = GetLocMatrix(LamToQ2key);
1480  }
1481 
1482  // Set up edge segment expansions from local geom info
1483  for (i = 0; i < nedges; ++i)
1484  {
1485  EdgeExp[i] = GetTraceExp(i);
1486  }
1487 
1488  // Set up matrix derived from <mu, Q_lam.n - \tau (U_lam - Lam) >
1489  for (i = 0; i < nbndry; ++i)
1490  {
1491  cnt = 0;
1492 
1493  Vmath::Zero(nbndry, lam, 1);
1494  lam[i] = 1.0;
1495  SetTraceToGeomOrientation(EdgeExp, lam);
1496 
1497  for (e = 0; e < nedges; ++e)
1498  {
1499  order_e = EdgeExp[e]->GetNcoeffs();
1500  nquad_e = EdgeExp[e]->GetNumPoints(0);
1501 
1502  normals = GetTraceNormal(e);
1503  edgedir = GetTraceOrient(e);
1504 
1505  work = Array<OneD, NekDouble>(nquad_e);
1506  varcoeff_work = Array<OneD, NekDouble>(nquad_e);
1507 
1508  GetTraceToElementMap(e, emap, sign, edgedir);
1509 
1510  StdRegions::VarCoeffType VarCoeff[3] = {
1513 
1514  // Q0 * n0 (BQ_0 terms)
1515  Array<OneD, NekDouble> edgeCoeffs(order_e);
1516  Array<OneD, NekDouble> edgePhys(nquad_e);
1517  for (j = 0; j < order_e; ++j)
1518  {
1519  edgeCoeffs[j] = sign[j] * (*LamToQ[0])(emap[j], i);
1520  }
1521 
1522  EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
1523  // @TODO Var coeffs
1524  // Multiply by variable coefficient
1525  // if ((x =
1526  // varcoeffs.find(VarCoeff[0]))
1527  // != varcoeffs.end())
1528  // {
1529  // GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
1530  // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
1531  // }
1532  if (mmf)
1533  {
1535  0, e, EdgeExp[e], normals, varcoeffs);
1536  Vmath::Vmul(nquad_e, ncdotMF, 1, edgePhys, 1, work, 1);
1537  }
1538  else
1539  {
1540  Vmath::Vmul(nquad_e, normals[0], 1, edgePhys, 1, work,
1541  1);
1542  }
1543 
1544  // Q1 * n1 (BQ_1 terms)
1545  for (j = 0; j < order_e; ++j)
1546  {
1547  edgeCoeffs[j] = sign[j] * (*LamToQ[1])(emap[j], i);
1548  }
1549 
1550  EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
1551 
1552  // @TODO var coeffs
1553  // Multiply by variable coefficients
1554  // if ((x =
1555  // varcoeffs.find(VarCoeff[1]))
1556  // != varcoeffs.end())
1557  // {
1558  // GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
1559  // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
1560  // }
1561 
1562  if (mmf)
1563  {
1565  1, e, EdgeExp[e], normals, varcoeffs);
1566  Vmath::Vvtvp(nquad_e, ncdotMF, 1, edgePhys, 1, work, 1,
1567  work, 1);
1568  }
1569  else
1570  {
1571  Vmath::Vvtvp(nquad_e, normals[1], 1, edgePhys, 1, work,
1572  1, work, 1);
1573  }
1574 
1575  // Q2 * n2 (BQ_2 terms)
1576  if (coordim == 3)
1577  {
1578  for (j = 0; j < order_e; ++j)
1579  {
1580  edgeCoeffs[j] = sign[j] * (*LamToQ[2])(emap[j], i);
1581  }
1582 
1583  EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
1584  // @TODO var coeffs
1585  // Multiply by variable coefficients
1586  // if ((x =
1587  // varcoeffs.find(VarCoeff[2]))
1588  // != varcoeffs.end())
1589  // {
1590  // GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
1591  // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
1592  // }
1593 
1594  Vmath::Vvtvp(nquad_e, normals[2], 1, edgePhys, 1, work,
1595  1, work, 1);
1596  }
1597 
1598  // - tau (ulam - lam)
1599  // Corresponds to the G and BU terms.
1600  for (j = 0; j < order_e; ++j)
1601  {
1602  edgeCoeffs[j] =
1603  sign[j] * LamToU(emap[j], i) - lam[cnt + j];
1604  }
1605 
1606  EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
1607 
1608  // Multiply by variable coefficients
1609  if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
1610  {
1612  e, EdgeExp[e], x->second, varcoeff_work);
1613  Vmath::Vmul(nquad_e, varcoeff_work, 1, edgePhys, 1,
1614  edgePhys, 1);
1615  }
1616 
1617  Vmath::Svtvp(nquad_e, -tau, edgePhys, 1, work, 1, work, 1);
1618  /// TODO: Add variable coeffs
1619  EdgeExp[e]->IProductWRTBase(work, edgeCoeffs);
1620 
1621  EdgeExp[e]->SetCoeffsToOrientation(edgedir, edgeCoeffs,
1622  edgeCoeffs);
1623 
1624  for (j = 0; j < order_e; ++j)
1625  {
1626  BndMat(cnt + j, i) = edgeCoeffs[j];
1627  }
1628 
1629  cnt += order_e;
1630  }
1631  }
1632  }
1633  break;
1634  // HDG postprocessing
1636  {
1638  mkey.GetConstFactors(), mkey.GetVarCoeffs());
1639  DNekScalMat &LapMat = *GetLocMatrix(lapkey);
1640 
1642  LapMat.GetRows(), LapMat.GetColumns());
1643  DNekMatSharedPtr lmat = returnval;
1644 
1645  (*lmat) = LapMat;
1646 
1647  // replace first column with inner product wrt 1
1648  int nq = GetTotPoints();
1649  Array<OneD, NekDouble> tmp(nq);
1651  Vmath::Fill(nq, 1.0, tmp, 1);
1652  IProductWRTBase(tmp, outarray);
1653 
1654  Vmath::Vcopy(m_ncoeffs, &outarray[0], 1, &(lmat->GetPtr())[0], 1);
1655  lmat->Invert();
1656  }
1657  break;
1659  {
1660  int ntraces = GetNtraces();
1661  int ncoords = GetCoordim();
1662  int nphys = GetTotPoints();
1664  Array<OneD, NekDouble> phys(nphys);
1665  returnval =
1667  DNekMat &Mat = *returnval;
1668  Vmath::Zero(m_ncoeffs * m_ncoeffs, Mat.GetPtr(), 1);
1669 
1671 
1672  for (int d = 0; d < ncoords; ++d)
1673  {
1674  Deriv[d] = Array<OneD, NekDouble>(nphys);
1675  }
1676 
1677  Array<OneD, int> tracepts(ntraces);
1678  Array<OneD, ExpansionSharedPtr> traceExp(ntraces);
1679  int maxtpts = 0;
1680  for (int t = 0; t < ntraces; ++t)
1681  {
1682  traceExp[t] = GetTraceExp(t);
1683  tracepts[t] = traceExp[t]->GetTotPoints();
1684  maxtpts = (maxtpts > tracepts[t]) ? maxtpts : tracepts[t];
1685  }
1686 
1687  Array<OneD, NekDouble> val(maxtpts), tmp, tmp1;
1688 
1689  Array<OneD, Array<OneD, NekDouble>> dphidn(ntraces);
1690  for (int t = 0; t < ntraces; ++t)
1691  {
1692  dphidn[t] =
1693  Array<OneD, NekDouble>(m_ncoeffs * tracepts[t], 0.0);
1694  }
1695 
1696  for (int i = 0; i < m_ncoeffs; ++i)
1697  {
1698  FillMode(i, phys);
1699  PhysDeriv(phys, Deriv[0], Deriv[1], Deriv[2]);
1700 
1701  for (int t = 0; t < ntraces; ++t)
1702  {
1703  const NormalVector norm = GetTraceNormal(t);
1704 
1706  LibUtilities::BasisKey tokey =
1707  traceExp[t]->GetBasis(0)->GetBasisKey();
1708  bool DoInterp = (fromkey != tokey);
1709 
1710  Array<OneD, NekDouble> n(tracepts[t]);
1711  ;
1712  for (int d = 0; d < ncoords; ++d)
1713  {
1714  // if variable p may need to interpolate
1715  if (DoInterp)
1716  {
1717  LibUtilities::Interp1D(fromkey, norm[d], tokey, n);
1718  }
1719  else
1720  {
1721  n = norm[d];
1722  }
1723 
1724  GetTracePhysVals(t, traceExp[t], Deriv[d], val,
1725  v_GetTraceOrient(t));
1726 
1727  Vmath::Vvtvp(tracepts[t], n, 1, val, 1,
1728  tmp = dphidn[t] + i * tracepts[t], 1,
1729  tmp1 = dphidn[t] + i * tracepts[t], 1);
1730  }
1731  }
1732  }
1733 
1734  for (int t = 0; t < ntraces; ++t)
1735  {
1736  int nt = tracepts[t];
1737  NekDouble h, p;
1738  TraceNormLen(t, h, p);
1739 
1740  // scaling from GJP paper
1741  NekDouble scale =
1742  (p == 1) ? 0.02 * h * h : 0.8 * pow(p + 1, -4.0) * h * h;
1743 
1744  for (int i = 0; i < m_ncoeffs; ++i)
1745  {
1746  for (int j = i; j < m_ncoeffs; ++j)
1747  {
1748  Vmath::Vmul(nt, dphidn[t] + i * nt, 1,
1749  dphidn[t] + j * nt, 1, val, 1);
1750  Mat(i, j) =
1751  Mat(i, j) + scale * traceExp[t]->Integral(val);
1752  }
1753  }
1754  }
1755 
1756  // fill in symmetric components.
1757  for (int i = 0; i < m_ncoeffs; ++i)
1758  {
1759  for (int j = 0; j < i; ++j)
1760  {
1761  Mat(i, j) = Mat(j, i);
1762  }
1763  }
1764  }
1765  break;
1766  default:
1767  ASSERTL0(false,
1768  "This matrix type cannot be generated from this class");
1769  break;
1770  }
1771 
1772  return returnval;
1773 }
1774 
1775 // Evaluate Coefficients of weak deriviative in the direction dir
1776 // given the input coefficicents incoeffs and the imposed
1777 // boundary values in EdgeExp (which will have its phys space updated);
1779  const Array<OneD, const NekDouble> &incoeffs,
1781  Array<OneD, Array<OneD, NekDouble>> &edgeCoeffs,
1782  Array<OneD, NekDouble> &out_d)
1783 {
1787 
1788  int ncoeffs = GetNcoeffs();
1789 
1791  DNekScalMat &Dmat = *GetLocMatrix(DerivType[dir]);
1792 
1793  Array<OneD, NekDouble> coeffs = incoeffs;
1794  DNekVec Coeffs(ncoeffs, coeffs, eWrapper);
1795 
1796  Coeffs = Transpose(Dmat) * Coeffs;
1797  Vmath::Neg(ncoeffs, coeffs, 1);
1798 
1799  // Add the boundary integral including the relevant part of
1800  // the normal
1801  AddNormTraceInt(dir, EdgeExp, edgeCoeffs, coeffs);
1802 
1803  DNekVec Out_d(ncoeffs, out_d, eWrapper);
1804 
1805  Out_d = InvMass * Coeffs;
1806 }
1807 
1809 {
1813 };
1814 
1816  const int edge, const Array<OneD, const NekDouble> &primCoeffs,
1817  DNekMatSharedPtr &inoutmat)
1818 {
1820  "Not set up for non boundary-interior expansions");
1821  ASSERTL1(inoutmat->GetRows() == inoutmat->GetColumns(),
1822  "Assuming that input matrix was square");
1823  int i, j;
1824  int id1, id2;
1825  ExpansionSharedPtr edgeExp = m_traceExp[edge].lock();
1826  int order_e = edgeExp->GetNcoeffs();
1827 
1830 
1831  StdRegions::VarCoeffMap varcoeffs;
1832  varcoeffs[StdRegions::eVarCoeffMass] = primCoeffs;
1833 
1836  varcoeffs);
1837  DNekScalMat &edgemat = *edgeExp->GetLocMatrix(mkey);
1838 
1839  // Now need to identify a map which takes the local edge
1840  // mass matrix to the matrix stored in inoutmat;
1841  // This can currently be deduced from the size of the matrix
1842 
1843  // - if inoutmat.m_rows() == v_NCoeffs() it is a full
1844  // matrix system
1845 
1846  // - if inoutmat.m_rows() == v_NumBndCoeffs() it is a
1847  // boundary CG system
1848 
1849  // - if inoutmat.m_rows() == v_NumDGBndCoeffs() it is a
1850  // trace DG system
1851  int rows = inoutmat->GetRows();
1852 
1853  if (rows == GetNcoeffs())
1854  {
1855  GetTraceToElementMap(edge, map, sign, v_GetTraceOrient(edge));
1856  }
1857  else if (rows == NumBndryCoeffs())
1858  {
1859  int nbndry = NumBndryCoeffs();
1860  Array<OneD, unsigned int> bmap(nbndry);
1861 
1862  GetTraceToElementMap(edge, map, sign, v_GetTraceOrient(edge));
1863 
1864  GetBoundaryMap(bmap);
1865 
1866  for (i = 0; i < order_e; ++i)
1867  {
1868  for (j = 0; j < nbndry; ++j)
1869  {
1870  if (map[i] == bmap[j])
1871  {
1872  map[i] = j;
1873  break;
1874  }
1875  }
1876  ASSERTL1(j != nbndry, "Did not find number in map");
1877  }
1878  }
1879  else if (rows == NumDGBndryCoeffs())
1880  {
1881  // possibly this should be a separate method
1882  int cnt = 0;
1883  map = Array<OneD, unsigned int>(order_e);
1884  sign = Array<OneD, int>(order_e, 1);
1885 
1886  for (i = 0; i < edge; ++i)
1887  {
1888  cnt += GetTraceNcoeffs(i);
1889  }
1890 
1891  for (i = 0; i < order_e; ++i)
1892  {
1893  map[i] = cnt++;
1894  }
1895  // check for mapping reversal
1897  {
1898  switch (edgeExp->GetBasis(0)->GetBasisType())
1899  {
1901  reverse(map.get(), map.get() + order_e);
1902  break;
1904  reverse(map.get(), map.get() + order_e);
1905  break;
1907  {
1908  swap(map[0], map[1]);
1909  for (i = 3; i < order_e; i += 2)
1910  {
1911  sign[i] = -1;
1912  }
1913  }
1914  break;
1915  default:
1916  ASSERTL0(false,
1917  "Edge boundary type not valid for this method");
1918  }
1919  }
1920  }
1921  else
1922  {
1923  ASSERTL0(false, "Could not identify matrix type from dimension");
1924  }
1925 
1926  for (i = 0; i < order_e; ++i)
1927  {
1928  id1 = map[i];
1929  for (j = 0; j < order_e; ++j)
1930  {
1931  id2 = map[j];
1932  (*inoutmat)(id1, id2) += edgemat(i, j) * sign[i] * sign[j];
1933  }
1934  }
1935 }
1936 
1937 /**
1938  * Given an edge and vector of element coefficients:
1939  * - maps those elemental coefficients corresponding to the edge into
1940  * an edge-vector.
1941  * - resets the element coefficients
1942  * - multiplies the edge vector by the edge mass matrix
1943  * - maps the edge coefficients back onto the elemental coefficients
1944  */
1946  const int edgeid, const Array<OneD, const NekDouble> &primCoeffs,
1947  const Array<OneD, NekDouble> &incoeffs, Array<OneD, NekDouble> &coeffs)
1948 {
1950  "Not set up for non boundary-interior expansions");
1951  int i;
1952  ExpansionSharedPtr edgeExp = m_traceExp[edgeid].lock();
1953  int order_e = edgeExp->GetNcoeffs();
1954 
1957 
1958  StdRegions::VarCoeffMap varcoeffs;
1959  varcoeffs[StdRegions::eVarCoeffMass] = primCoeffs;
1960 
1963  varcoeffs);
1964  DNekScalMat &edgemat = *edgeExp->GetLocMatrix(mkey);
1965 
1966  NekVector<NekDouble> vEdgeCoeffs(order_e);
1967 
1968  GetTraceToElementMap(edgeid, map, sign, v_GetTraceOrient(edgeid));
1969 
1970  for (i = 0; i < order_e; ++i)
1971  {
1972  vEdgeCoeffs[i] = incoeffs[map[i]] * sign[i];
1973  }
1974 
1975  vEdgeCoeffs = edgemat * vEdgeCoeffs;
1976 
1977  for (i = 0; i < order_e; ++i)
1978  {
1979  coeffs[map[i]] += vEdgeCoeffs[i] * sign[i];
1980  }
1981 }
1982 
1984  const DNekScalMatSharedPtr &r_bnd)
1985 {
1986  MatrixStorage storage = eFULL;
1987  DNekMatSharedPtr m_vertexmatrix;
1988 
1989  int nVerts, vid1, vid2, vMap1, vMap2;
1990  NekDouble VertexValue;
1991 
1992  nVerts = GetNverts();
1993 
1994  m_vertexmatrix =
1995  MemoryManager<DNekMat>::AllocateSharedPtr(nVerts, nVerts, 0.0, storage);
1996  DNekMat &VertexMat = (*m_vertexmatrix);
1997 
1998  for (vid1 = 0; vid1 < nVerts; ++vid1)
1999  {
2000  vMap1 = GetVertexMap(vid1);
2001 
2002  for (vid2 = 0; vid2 < nVerts; ++vid2)
2003  {
2004  vMap2 = GetVertexMap(vid2);
2005  VertexValue = (*r_bnd)(vMap1, vMap2);
2006  VertexMat.SetValue(vid1, vid2, VertexValue);
2007  }
2008  }
2009 
2010  return m_vertexmatrix;
2011 }
2012 
2013 void Expansion2D::v_GenTraceExp(const int traceid, ExpansionSharedPtr &exp)
2014 {
2016  GetTraceBasisKey(traceid), m_geom->GetEdge(traceid));
2017 }
2018 
2020 {
2021  int n, j;
2022  int nEdgeCoeffs;
2023  int nBndCoeffs = NumBndryCoeffs();
2024 
2025  Array<OneD, unsigned int> bmap(nBndCoeffs);
2026  GetBoundaryMap(bmap);
2027 
2028  // Map from full system to statically condensed system (i.e reverse
2029  // GetBoundaryMap)
2030  map<int, int> invmap;
2031  for (j = 0; j < nBndCoeffs; ++j)
2032  {
2033  invmap[bmap[j]] = j;
2034  }
2035 
2036  // Number of interior edge coefficients
2037  nEdgeCoeffs = GetTraceNcoeffs(eid) - 2;
2038 
2040 
2041  Array<OneD, unsigned int> edgemaparray(nEdgeCoeffs);
2042  Array<OneD, unsigned int> maparray(nEdgeCoeffs);
2043  Array<OneD, int> signarray(nEdgeCoeffs, 1);
2044  StdRegions::Orientation eOrient = geom->GetEorient(eid);
2045 
2046  // maparray is the location of the edge within the matrix
2047  GetTraceInteriorToElementMap(eid, maparray, signarray, eOrient);
2048 
2049  for (n = 0; n < nEdgeCoeffs; ++n)
2050  {
2051  edgemaparray[n] = invmap[maparray[n]];
2052  }
2053 
2054  return edgemaparray;
2055 }
2056 
2058 {
2059  v_ComputeTraceNormal(edge);
2060 }
2061 
2063  Array<OneD, int> &idmap, const int nq0,
2064  const int nq1)
2065 {
2066  boost::ignore_unused(nq1);
2067 
2068  if (idmap.size() != nq0)
2069  {
2070  idmap = Array<OneD, int>(nq0);
2071  }
2072  switch (orient)
2073  {
2074  case StdRegions::eForwards:
2075  // Fwd
2076  for (int i = 0; i < nq0; ++i)
2077  {
2078  idmap[i] = i;
2079  }
2080  break;
2082  {
2083  // Bwd
2084  for (int i = 0; i < nq0; ++i)
2085  {
2086  idmap[i] = nq0 - 1 - i;
2087  }
2088  }
2089  break;
2090  default:
2091  ASSERTL0(false, "Unknown orientation");
2092  break;
2093  }
2094 }
2095 
2097  const int dir, const int shapedim, const StdRegions::VarCoeffMap &varcoeffs)
2098 {
2099  return Expansion::v_GetMF(dir, shapedim, varcoeffs);
2100 }
2101 
2103  const int dir, const StdRegions::VarCoeffMap &varcoeffs)
2104 {
2105  return Expansion::v_GetMFDiv(dir, varcoeffs);
2106 }
2107 
2109  const int dir, const StdRegions::VarCoeffMap &varcoeffs)
2110 {
2111  return Expansion::v_GetMFMag(dir, varcoeffs);
2112 }
2113 
2114 // Compute edgenormal \cdot vector
2116  const int dir, const int edge, ExpansionSharedPtr &EdgeExp_e,
2117  const Array<OneD, const Array<OneD, NekDouble>> &normals,
2118  const StdRegions::VarCoeffMap &varcoeffs)
2119 {
2120  int nquad_e = EdgeExp_e->GetNumPoints(0);
2121  int coordim = GetCoordim();
2122  int nquad0 = m_base[0]->GetNumPoints();
2123  int nquad1 = m_base[1]->GetNumPoints();
2124  int nqtot = nquad0 * nquad1;
2125 
2126  StdRegions::VarCoeffType MMFCoeffs[15] = {
2135 
2136  StdRegions::VarCoeffMap::const_iterator MFdir;
2137 
2138  Array<OneD, NekDouble> ncdotMF(nqtot, 0.0);
2139  Array<OneD, NekDouble> tmp(nqtot);
2140  Array<OneD, NekDouble> tmp_e(nquad_e);
2141  for (int k = 0; k < coordim; k++)
2142  {
2143  MFdir = varcoeffs.find(MMFCoeffs[dir * 5 + k]);
2144  tmp = MFdir->second;
2145 
2146  GetPhysEdgeVarCoeffsFromElement(edge, EdgeExp_e, tmp, tmp_e);
2147 
2148  Vmath::Vvtvp(nquad_e, &tmp_e[0], 1, &normals[k][0], 1, &ncdotMF[0], 1,
2149  &ncdotMF[0], 1);
2150  }
2151  return ncdotMF;
2152 }
2153 
2155  const Array<OneD, Array<OneD, NekDouble>> &vec)
2156 {
2158  GetLeftAdjacentElementExp()->GetTraceNormal(
2160 
2161  int nq = GetTotPoints();
2162  Array<OneD, NekDouble> Fn(nq);
2163  Vmath::Vmul(nq, &vec[0][0], 1, &normals[0][0], 1, &Fn[0], 1);
2164  Vmath::Vvtvp(nq, &vec[1][0], 1, &normals[1][0], 1, &Fn[0], 1, &Fn[0], 1);
2165  Vmath::Vvtvp(nq, &vec[2][0], 1, &normals[2][0], 1, &Fn[0], 1, &Fn[0], 1);
2166 
2167  return StdExpansion::Integral(Fn);
2168 }
2169 
2170 void Expansion2D::v_TraceNormLen(const int traceid, NekDouble &h, NekDouble &p)
2171 {
2173 
2174  int nverts = geom->GetNumVerts();
2175 
2176  // vertices on edges
2177  SpatialDomains::PointGeom ev0 = *geom->GetVertex(traceid);
2178  SpatialDomains::PointGeom ev1 = *geom->GetVertex((traceid + 1) % nverts);
2179 
2180  // vertex on adjacent edge to ev0
2182  *geom->GetVertex((traceid + (nverts - 1)) % nverts);
2183 
2184  // calculate perpendicular distance of normal length
2185  // from first vertex
2186  NekDouble h1 = ev0.dist(vadj);
2187  SpatialDomains::PointGeom Dx, Dx1;
2188 
2189  Dx.Sub(ev1, ev0);
2190  Dx1.Sub(vadj, ev0);
2191 
2192  NekDouble d1 = Dx.dot(Dx1);
2193  NekDouble lenDx = Dx.dot(Dx);
2194  h = sqrt(h1 * h1 - d1 * d1 / lenDx);
2195 
2196  // perpendicular distanace from second vertex
2197  SpatialDomains::PointGeom vadj1 = *geom->GetVertex((traceid + 2) % nverts);
2198 
2199  h1 = ev1.dist(vadj1);
2200  Dx1.Sub(vadj1, ev1);
2201  d1 = Dx.dot(Dx1);
2202 
2203  h = (h + sqrt(h1 * h1 - d1 * d1 / lenDx)) * 0.5;
2204 
2205  int dirn = (geom->GetDir(traceid) == 0) ? 1 : 0;
2206 
2207  p = (NekDouble)(GetBasisNumModes(dirn) - 1);
2208 }
2209 } // namespace LocalRegions
2210 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:272
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:15
Describes the specification for a Basis.
Definition: Basis.h:50
virtual void v_AddRobinTraceContribution(const int traceid, const Array< OneD, const NekDouble > &primCoeffs, const Array< OneD, NekDouble > &incoeffs, Array< OneD, NekDouble > &coeffs)
std::vector< bool > m_requireNeg
Definition: Expansion2D.h:119
Array< OneD, NekDouble > v_GetnEdgecdotMF(const int dir, const int edge, ExpansionSharedPtr &EdgeExp_e, const Array< OneD, const Array< OneD, NekDouble >> &normals, const StdRegions::VarCoeffMap &varcoeffs)
virtual Array< OneD, NekDouble > v_GetMFDiv(const int dir, const StdRegions::VarCoeffMap &varcoeffs)
void SetTraceToGeomOrientation(Array< OneD, ExpansionSharedPtr > &EdgeExp, Array< OneD, NekDouble > &inout)
virtual void v_GenTraceExp(const int traceid, ExpansionSharedPtr &exp)
void v_ReOrientTracePhysMap(const StdRegions::Orientation orient, Array< OneD, int > &idmap, const int nq0, const int nq1)
virtual void v_AddEdgeNormBoundaryInt(const int edge, const ExpansionSharedPtr &EdgeExp, const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
DNekScalMatSharedPtr CreateMatrix(const MatrixKey &mkey)
Definition: Expansion2D.cpp:59
virtual Array< OneD, NekDouble > v_GetMFMag(const int dir, const StdRegions::VarCoeffMap &varcoeffs)
Array< OneD, unsigned int > GetTraceInverseBoundaryMap(int eid)
virtual void v_SetUpPhysNormals(const int edge)
virtual void v_TraceNormLen(const int traceid, NekDouble &h, NekDouble &p)
virtual void v_AddRobinMassMatrix(const int edgeid, const Array< OneD, const NekDouble > &primCoeffs, DNekMatSharedPtr &inoutmat)
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
virtual Array< OneD, NekDouble > v_GetMF(const int dir, const int shapedim, const StdRegions::VarCoeffMap &varcoeffs)
SpatialDomains::Geometry2DSharedPtr GetGeom2D() const
Definition: Expansion2D.h:180
void AddHDGHelmholtzEdgeTerms(const NekDouble tau, const int edge, Array< OneD, ExpansionSharedPtr > &EdgeExp, Array< OneD, NekDouble > &edgePhys, const StdRegions::VarCoeffMap &dirForcing, Array< OneD, NekDouble > &outarray)
void AddNormTraceInt(const int dir, Array< OneD, ExpansionSharedPtr > &EdgeExp, Array< OneD, Array< OneD, NekDouble >> &edgeCoeffs, Array< OneD, NekDouble > &outarray)
void GetPhysEdgeVarCoeffsFromElement(const int edge, ExpansionSharedPtr &EdgeExp, const Array< OneD, const NekDouble > &varcoeff, Array< OneD, NekDouble > &outarray)
void AddEdgeBoundaryInt(const int edge, ExpansionSharedPtr &EdgeExp, Array< OneD, NekDouble > &edgePhys, Array< OneD, NekDouble > &outarray, const StdRegions::VarCoeffMap &varcoeffs=StdRegions::NullVarCoeffMap)
virtual DNekMatSharedPtr v_BuildVertexMatrix(const DNekScalMatSharedPtr &r_bnd)
virtual NekDouble v_VectorFlux(const Array< OneD, Array< OneD, NekDouble >> &vec)
virtual void v_DGDeriv(const int dir, const Array< OneD, const NekDouble > &incoeffs, Array< OneD, ExpansionSharedPtr > &EdgeExp, Array< OneD, Array< OneD, NekDouble >> &edgeCoeffs, Array< OneD, NekDouble > &out_d)
void AddHDGHelmholtzTraceTerms(const NekDouble tau, const Array< OneD, const NekDouble > &inarray, Array< OneD, ExpansionSharedPtr > &EdgeExp, const StdRegions::VarCoeffMap &dirForcing, Array< OneD, NekDouble > &outarray)
SpatialDomains::GeometrySharedPtr GetGeom() const
Definition: Expansion.cpp:166
SpatialDomains::GeometrySharedPtr m_geom
Definition: Expansion.h:272
Array< OneD, NekDouble > v_GetMF(const int dir, const int shapedim, const StdRegions::VarCoeffMap &varcoeffs)
Definition: Expansion.cpp:627
ExpansionSharedPtr GetLeftAdjacentElementExp() const
Definition: Expansion.h:433
DNekMatSharedPtr BuildVertexMatrix(const DNekScalMatSharedPtr &r_bnd)
Definition: Expansion.cpp:100
void ComputeGmatcdotMF(const Array< TwoD, const NekDouble > &df, const Array< OneD, const NekDouble > &direction, Array< OneD, Array< OneD, NekDouble >> &dfdir)
Definition: Expansion.cpp:597
virtual void v_ComputeTraceNormal(const int id)
Definition: Expansion.cpp:861
int GetLeftAdjacentElementTrace() const
Definition: Expansion.h:446
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:199
Array< OneD, NekDouble > v_GetMFMag(const int dir, const StdRegions::VarCoeffMap &varcoeffs)
Definition: Expansion.cpp:703
std::map< int, ExpansionWeakPtr > m_traceExp
Definition: Expansion.h:271
Array< OneD, NekDouble > v_GetMFDiv(const int dir, const StdRegions::VarCoeffMap &varcoeffs)
Definition: Expansion.cpp:680
void AddEdgeNormBoundaryInt(const int edge, const std::shared_ptr< Expansion > &EdgeExp, const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
Definition: Expansion.cpp:113
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
Definition: Expansion.h:273
virtual StdRegions::Orientation v_GetTraceOrient(int trace)
Definition: Expansion.cpp:806
ExpansionSharedPtr GetTraceExp(const int traceid)
Definition: Expansion.h:406
StdRegions::Orientation GetTraceOrient(int trace)
Definition: Expansion.h:167
IndexMapValuesSharedPtr GetIndexMap(const IndexMapKey &ikey)
Definition: Expansion.h:145
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition: Expansion.cpp:88
void TraceNormLen(const int traceid, NekDouble &h, NekDouble &p)
Definition: Expansion.h:249
const NormalVector & GetTraceNormal(const int id)
Definition: Expansion.cpp:250
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
PointGeomSharedPtr GetVertex(int i) const
Returns vertex i of this object.
Definition: Geometry.h:345
void Sub(PointGeom &a, PointGeom &b)
Definition: PointGeom.cpp:124
NekDouble dot(PointGeom &a)
retun the dot product between this and input a
Definition: PointGeom.cpp:187
NekDouble dist(PointGeom &a)
return distance between this and input a
Definition: PointGeom.cpp:180
void GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:677
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:130
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:140
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:499
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:611
const LibUtilities::PointsKeyVector GetPointsKeys() const
DNekScalBlkMatSharedPtr GetLocStaticCondMatrix(const LocalRegions::MatrixKey &mkey)
Definition: StdExpansion.h:649
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:307
int GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
Definition: StdExpansion.h:687
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:536
int GetNtraces() const
Returns the number of trace elements connected to this element.
Definition: StdExpansion.h:359
int GetNverts() const
This function returns the number of vertices of the expansion domain.
Definition: StdExpansion.h:254
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:692
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:375
void GetTraceInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eForwards)
Definition: StdExpansion.h:716
int GetTraceNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th trace.
Definition: StdExpansion.h:269
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:845
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space.
LibUtilities::NekManager< StdMatrixKey, DNekMat, StdMatrixKey::opLess > m_stdMatrixManager
int GetBasisNumModes(const int dir) const
This function returns the number of expansion modes in the dir direction.
Definition: StdExpansion.h:176
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:850
Array< OneD, LibUtilities::BasisSharedPtr > m_base
const ConstFactorMap & GetConstFactors() const
Definition: StdMatrixKey.h:140
LibUtilities::ShapeType GetShapeType() const
Definition: StdMatrixKey.h:90
const VarCoeffMap & GetVarCoeffs() const
Definition: StdMatrixKey.h:167
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:85
const Array< OneD, const NekDouble > & GetVarCoeff(const StdRegions::VarCoeffType &coeff) const
Definition: StdMatrixKey.h:150
bool HasVarCoeff(const StdRegions::VarCoeffType &coeff) const
Definition: StdMatrixKey.h:172
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:126
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:135
const VarCoeffMap GetVarCoeffAsMap(const VarCoeffType &coeff) const
Definition: StdMatrixKey.h:160
void Interp1D(const BasisKey &fbasis0, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, Array< OneD, NekDouble > &to)
this function interpolates a 1D function evaluated at the quadrature points of the basis fbasis0 to ...
Definition: Interp.cpp:52
void InterpCoeff1D(const BasisKey &fbasis0, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, Array< OneD, NekDouble > &to)
Definition: InterpCoeff.cpp:46
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:250
@ eGauss_Lagrange
Lagrange Polynomials using the Gauss points.
Definition: BasisType.h:59
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:58
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:50
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
std::shared_ptr< IndexMapValues > IndexMapValuesSharedPtr
Definition: IndexMapKey.h:128
std::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:51
@ eNoGeomType
No type defined.
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry.h:65
std::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:53
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
Definition: StdRegions.hpp:240
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:282
static ConstFactorMap NullConstFactorMap
Definition: StdRegions.hpp:283
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:79
static Array< OneD, NekDouble > NullNekDouble1DArray
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
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:209
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:622
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
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:574
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291