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 
37 #include <LocalRegions/Expansion2D.h>
38 #include <LocalRegions/Expansion1D.h>
39 #include <LocalRegions/Expansion3D.h>
40 #include <SpatialDomains/Geometry.h>
41 #include <SpatialDomains/Geometry2D.h>
42 #include <LibUtilities/Foundations/InterpCoeff.h>
43 #include <LocalRegions/MatrixKey.h>
44 #include <LibUtilities/Foundations/ManagerAccess.h>
45 
46 using namespace std;
47 
48 namespace Nektar
49 {
50  namespace LocalRegions
51  {
52  Expansion2D::Expansion2D(SpatialDomains::Geometry2DSharedPtr pGeom) :
53  StdExpansion(), Expansion(pGeom), StdExpansion2D()
54  {
55  m_elementFaceLeft = -1;
56  m_elementFaceRight = -1;
57  }
58 
59  void Expansion2D::v_AddEdgeNormBoundaryInt(
60  const int edge,
61  const ExpansionSharedPtr &EdgeExp,
62  const Array<OneD, const NekDouble> &Fx,
63  const Array<OneD, const NekDouble> &Fy,
64  Array<OneD, NekDouble> &outarray)
65  {
66  ASSERTL1(GetCoordim() == 2,
67  "Routine only set up for two-dimensions");
68 
69  const Array<OneD, const Array<OneD, NekDouble> > normals
70  = GetEdgeNormal(edge);
71 
72  if (m_requireNeg.size() == 0)
73  {
74  m_requireNeg.resize(GetNedges());
75 
76  for (int i = 0; i < GetNedges(); ++i)
77  {
78  m_requireNeg[i] = false;
79  if (m_negatedNormals[i])
80  {
81  m_requireNeg[i] = true;
82  continue;
83  }
84 
85  Expansion1DSharedPtr edgeExp =
86  m_edgeExp[i].lock()->as<Expansion1D>();
87 
88  if (edgeExp->GetRightAdjacentElementExp())
89  {
90  if (edgeExp->GetRightAdjacentElementExp()->GetGeom2D()
91  ->GetGlobalID() == GetGeom2D()->GetGlobalID())
92  {
93  m_requireNeg[i] = true;
94  }
95  }
96  }
97  }
98 
99  // We allow the case of mixed polynomial order by supporting only
100  // those modes on the edge common to both adjoining elements. This
101  // is enforced here by taking the minimum size and padding with
102  // zeros.
103  int nquad_e = min(EdgeExp->GetNumPoints(0),
104  int(normals[0].num_elements()));
105 
106  int nEdgePts = EdgeExp->GetTotPoints();
107  Array<OneD, NekDouble> edgePhys(nEdgePts);
108  Vmath::Vmul (nquad_e, normals[0], 1, Fx, 1, edgePhys, 1);
109  Vmath::Vvtvp(nquad_e, normals[1], 1, Fy, 1, edgePhys, 1,
110  edgePhys, 1);
111 
112  Expansion1DSharedPtr locExp = EdgeExp->as<Expansion1D>();
113 
114  if (m_negatedNormals[edge])
115  {
116  Vmath::Neg(nquad_e, edgePhys, 1);
117  }
118  else if (locExp->GetRightAdjacentElementEdge() != -1)
119  {
120  if (locExp->GetRightAdjacentElementExp()->GetGeom2D()->GetGlobalID()
121  == GetGeom2D()->GetGlobalID())
122  {
123  Vmath::Neg(nquad_e, edgePhys, 1);
124  }
125  }
126 
127  AddEdgeNormBoundaryInt(edge, EdgeExp, edgePhys, outarray);
128  }
129 
130  void Expansion2D::v_AddEdgeNormBoundaryInt(
131  const int edge,
132  const ExpansionSharedPtr &EdgeExp,
133  const Array<OneD, const NekDouble> &Fn,
134  Array<OneD, NekDouble> &outarray)
135  {
136  int i;
137 
138  if (m_requireNeg.size() == 0)
139  {
140  m_requireNeg.resize(GetNedges());
141 
142  for (i = 0; i < GetNedges(); ++i)
143  {
144  m_requireNeg[i] = false;
145  if (m_negatedNormals[i])
146  {
147  m_requireNeg[i] = true;
148  continue;
149  }
150 
151  Expansion1DSharedPtr edgeExp =
152  m_edgeExp[i].lock()->as<Expansion1D>();
153 
154  if (edgeExp->GetRightAdjacentElementExp())
155  {
156  if (edgeExp->GetRightAdjacentElementExp()->GetGeom2D()
157  ->GetGlobalID() == GetGeom2D()->GetGlobalID())
158  {
159  m_requireNeg[i] = true;
160  }
161  }
162  }
163  }
164 
165  StdRegions::IndexMapKey ikey(
166  StdRegions::eEdgeToElement, DetShapeType(),
167  GetBasisNumModes(0), GetBasisNumModes(1), 0,
168  edge, GetEorient(edge));
169  StdRegions::IndexMapValuesSharedPtr map =
170  StdExpansion::GetIndexMap(ikey);
171 
172  // Order of the element
173  int order_e = map->num_elements();
174  // Order of the trace
175  int n_coeffs = EdgeExp->GetNcoeffs();
176 
177  Array<OneD, NekDouble> edgeCoeffs(n_coeffs);
178  if(n_coeffs!=order_e) // Going to orthogonal space
179  {
180  EdgeExp->FwdTrans(Fn, edgeCoeffs);
181  Expansion1DSharedPtr locExp = EdgeExp->as<Expansion1D>();
182 
183  if (m_requireNeg[edge])
184  {
185  Vmath::Neg(n_coeffs, edgeCoeffs, 1);
186  }
187 
188  Array<OneD, NekDouble> coeff(n_coeffs,0.0);
189  LibUtilities::BasisType btype = ((LibUtilities::BasisType) 1); //1-->Ortho_A
190  LibUtilities::BasisKey bkey_ortho(btype,EdgeExp->GetBasis(0)->GetNumModes(),EdgeExp->GetBasis(0)->GetPointsKey());
191  LibUtilities::BasisKey bkey(EdgeExp->GetBasis(0)->GetBasisType(),EdgeExp->GetBasis(0)->GetNumModes(),EdgeExp->GetBasis(0)->GetPointsKey());
192  LibUtilities::InterpCoeff1D(bkey,edgeCoeffs,bkey_ortho,coeff);
193 
194  // Cutting high frequencies
195  for(i = order_e; i < n_coeffs; i++)
196  {
197  coeff[i] = 0.0;
198  }
199 
200  LibUtilities::InterpCoeff1D(bkey_ortho,coeff,bkey,edgeCoeffs);
201 
202  StdRegions::StdMatrixKey masskey(StdRegions::eMass,LibUtilities::eSegment,*EdgeExp);
203  EdgeExp->MassMatrixOp(edgeCoeffs, edgeCoeffs, masskey);
204  }
205  else
206  {
207  EdgeExp->IProductWRTBase(Fn, edgeCoeffs);
208 
209  Expansion1DSharedPtr locExp = EdgeExp->as<Expansion1D>();
210 
211  if (m_requireNeg[edge])
212  {
213  Vmath::Neg(n_coeffs, edgeCoeffs, 1);
214  }
215  }
216 
217  // Implementation for all the basis except Gauss points
218  if(EdgeExp->GetBasis(0)->GetBasisType() !=
219  LibUtilities::eGauss_Lagrange)
220  {
221  // add data to outarray if forward edge normal is outwards
222  for(i = 0; i < order_e; ++i)
223  {
224  outarray[(*map)[i].index] +=
225  (*map)[i].sign * edgeCoeffs[i];
226  }
227  }
228  else
229  {
230  int nCoeffs0, nCoeffs1;
231  int j;
232 
233  StdRegions::ConstFactorMap factors;
234  factors[StdRegions::eFactorGaussEdge] = edge;
235  StdRegions::StdMatrixKey key(StdRegions::eGaussDG,
236  DetShapeType(),*this,factors);
237 
238  DNekMatSharedPtr mat_gauss = m_stdMatrixManager[key];
239 
240  switch(edge)
241  {
242  case 0:
243  {
244  nCoeffs1 = m_base[1]->GetNumModes();
245 
246  for(i = 0; i < order_e; ++i)
247  {
248  for(j = 0; j < nCoeffs1; j++)
249  {
250  outarray[(*map)[i].index + j*order_e] +=
251  mat_gauss->GetPtr()[j]*
252  (*map)[i].sign*edgeCoeffs[i];
253  }
254  }
255  break;
256  }
257  case 1:
258  {
259  nCoeffs0 = m_base[0]->GetNumModes();
260 
261  for(i = 0; i < order_e; ++i)
262  {
263  for(j = 0; j < nCoeffs0; j++)
264  {
265  outarray[(*map)[i].index - j] +=
266  mat_gauss->GetPtr()[order_e - 1 -j]*
267  (*map)[i].sign*edgeCoeffs[i];
268  }
269  }
270  break;
271  }
272  case 2:
273  {
274  nCoeffs1 = m_base[1]->GetNumModes();
275 
276  for(i = 0; i < order_e; ++i)
277  {
278  for(j = 0; j < nCoeffs1; j++)
279  {
280  outarray[(*map)[i].index - j*order_e] +=
281  mat_gauss->GetPtr()[order_e - 1 - j]*
282  (*map)[i].sign*edgeCoeffs[i];
283  }
284  }
285  break;
286  }
287  case 3:
288  {
289  nCoeffs0 = m_base[0]->GetNumModes();
290 
291  for(i = 0; i < order_e; ++i)
292  {
293  for(j = 0; j < nCoeffs0; j++)
294  {
295  outarray[(*map)[i].index + j] +=
296  mat_gauss->GetPtr()[j]*
297  (*map)[i].sign*edgeCoeffs[i];
298  }
299  }
300  break;
301  }
302  default:
303  ASSERTL0(false,"edge value (< 3) is out of range");
304  break;
305  }
306  }
307  }
308 
309  void Expansion2D::SetTraceToGeomOrientation(
310  Array<OneD, ExpansionSharedPtr> &EdgeExp,
311  Array<OneD, NekDouble> &inout)
312  {
313  int i, cnt = 0;
314  int nedges = GetNedges();
315  Array<OneD, NekDouble> e_tmp;
316 
317  for(i = 0; i < nedges; ++i)
318  {
319  EdgeExp[i]->SetCoeffsToOrientation(GetEorient(i),
320  e_tmp = inout + cnt,
321  e_tmp = inout + cnt);
322  cnt += GetEdgeNcoeffs(i);
323  }
324  }
325 
326  /**
327  * Computes the C matrix entries due to the presence of the identity
328  * matrix in Eqn. 32.
329  */
330  void Expansion2D::AddNormTraceInt(
331  const int dir,
332  Array<OneD, const NekDouble> &inarray,
333  Array<OneD, ExpansionSharedPtr> &EdgeExp,
334  Array<OneD, NekDouble> &outarray,
335  const StdRegions::VarCoeffMap &varcoeffs)
336  {
337  int i,e,cnt;
338  int order_e,nquad_e;
339  int nedges = GetNedges();
340 
341  cnt = 0;
342  for(e = 0; e < nedges; ++e)
343  {
344  order_e = EdgeExp[e]->GetNcoeffs();
345  nquad_e = EdgeExp[e]->GetNumPoints(0);
346 
347  const Array<OneD, const Array<OneD, NekDouble> > &normals
348  = GetEdgeNormal(e);
349  Array<OneD, NekDouble> edgeCoeffs(order_e);
350  Array<OneD, NekDouble> edgePhys (nquad_e);
351 
352  for(i = 0; i < order_e; ++i)
353  {
354  edgeCoeffs[i] = inarray[i+cnt];
355  }
356  cnt += order_e;
357 
358  EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
359 
360  // Multiply by variable coefficient
361  /// @TODO: Document this
362  // StdRegions::VarCoeffType VarCoeff[3] = {StdRegions::eVarCoeffD00,
363  // StdRegions::eVarCoeffD11,
364  // StdRegions::eVarCoeffD22};
365  // StdRegions::VarCoeffMap::const_iterator x;
366  // Array<OneD, NekDouble> varcoeff_work(nquad_e);
367 
368  // if ((x = varcoeffs.find(VarCoeff[dir])) != varcoeffs.end())
369  // {
370  // GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
371  // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
372  // }
373 
374  if (varcoeffs.find(StdRegions::eVarCoeffMF1x) !=
375  varcoeffs.end())
376  {
377  // MMF case
378  Array<OneD, NekDouble> ncdotMF_e =
379  v_GetnEdgecdotMF(dir, e, EdgeExp[e], normals,
380  varcoeffs);
381 
382  Vmath::Vmul(nquad_e, ncdotMF_e, 1,
383  edgePhys, 1,
384  edgePhys, 1);
385  }
386  else
387  {
388  Vmath::Vmul(nquad_e, normals[dir], 1,
389  edgePhys, 1,
390  edgePhys, 1);
391  }
392 
393  if (m_negatedNormals[e])
394  {
395  Vmath::Neg(nquad_e, edgePhys, 1);
396  }
397 
398  AddEdgeBoundaryInt(e, EdgeExp[e], edgePhys, outarray, varcoeffs);
399  }
400  }
401 
402  void Expansion2D::AddNormTraceInt(
403  const int dir,
404  Array<OneD, ExpansionSharedPtr> &EdgeExp,
405  Array<OneD, Array<OneD, NekDouble> > &edgeCoeffs,
406  Array<OneD, NekDouble> &outarray)
407  {
408  int e;
409  int nquad_e;
410  int nedges = GetNedges();
411 
412  for(e = 0; e < nedges; ++e)
413  {
414  nquad_e = EdgeExp[e]->GetNumPoints(0);
415 
416  Array<OneD, NekDouble> edgePhys(nquad_e);
417  const Array<OneD, const Array<OneD, NekDouble> > &normals
418  = GetEdgeNormal(e);
419 
420  EdgeExp[e]->BwdTrans(edgeCoeffs[e], edgePhys);
421 
422  Vmath::Vmul(nquad_e, normals[dir], 1, edgePhys, 1, edgePhys, 1);
423 
424  if (m_negatedNormals[e])
425  {
426  Vmath::Neg(nquad_e, edgePhys, 1);
427  }
428 
429  AddEdgeBoundaryInt(e, EdgeExp[e], edgePhys, outarray);
430  }
431  }
432 
433 
434  /**
435  * For a given edge add the \tilde{F}_1j contributions
436  */
437  void Expansion2D::AddEdgeBoundaryInt(
438  const int edge,
439  ExpansionSharedPtr &EdgeExp,
440  Array<OneD, NekDouble> &edgePhys,
441  Array<OneD, NekDouble> &outarray,
442  const StdRegions::VarCoeffMap &varcoeffs)
443  {
444  int i;
445  int order_e = EdgeExp->GetNcoeffs();
446  int nquad_e = EdgeExp->GetNumPoints(0);
447  Array<OneD,unsigned int> map;
448  Array<OneD,int> sign;
449  Array<OneD, NekDouble> coeff(order_e);
450 
451  GetEdgeToElementMap(edge, v_GetEorient(edge), map, sign);
452 
453  StdRegions::VarCoeffType VarCoeff[3] = {StdRegions::eVarCoeffD00,
454  StdRegions::eVarCoeffD11,
455  StdRegions::eVarCoeffD22};
456  StdRegions::VarCoeffMap::const_iterator x;
457 
458  /// @TODO Variable coeffs
459  if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
460  {
461  Array<OneD, NekDouble> work(nquad_e);
462  GetPhysEdgeVarCoeffsFromElement(
463  edge, EdgeExp, x->second, work);
464  Vmath::Vmul(nquad_e, work, 1, edgePhys, 1, edgePhys, 1);
465  }
466 
467  EdgeExp->IProductWRTBase(edgePhys, coeff);
468 
469  // add data to out array
470  for(i = 0; i < order_e; ++i)
471  {
472  outarray[map[i]] += sign[i]*coeff[i];
473  }
474  }
475 
476  // This method assumes that data in EdgeExp is orientated according to
477  // elemental counter clockwise format AddHDGHelmholtzTraceTerms with
478  // directions
479  void Expansion2D::AddHDGHelmholtzTraceTerms(
480  const NekDouble tau,
481  const Array<OneD, const NekDouble> &inarray,
482  Array<OneD, ExpansionSharedPtr> &EdgeExp,
483  const StdRegions::VarCoeffMap &dirForcing,
484  Array<OneD,NekDouble> &outarray)
485  {
486  ASSERTL0(&inarray[0] != &outarray[0],
487  "Input and output arrays use the same memory");
488 
489  int e, cnt, order_e, nedges = GetNedges();
490  Array<OneD, const NekDouble> tmp;
491 
492  cnt = 0;
493 
494  for(e = 0; e < nedges; ++e)
495  {
496  order_e = EdgeExp[e]->GetNcoeffs();
497  Array<OneD, NekDouble> edgeCoeffs(order_e);
498  Array<OneD, NekDouble> edgePhys (EdgeExp[e]->GetTotPoints());
499 
500  Vmath::Vcopy(order_e, tmp = inarray + cnt, 1, edgeCoeffs, 1);
501  EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
502  AddHDGHelmholtzEdgeTerms(
503  tau, e, EdgeExp, edgePhys, dirForcing, outarray);
504 
505  cnt += order_e;
506  }
507  }
508 
509  // evaluate additional terms in HDG edges. Not that this assumes that
510  // edges are unpacked into local cartesian order.
511  void Expansion2D::AddHDGHelmholtzEdgeTerms(
512  const NekDouble tau,
513  const int edge,
514  Array<OneD, ExpansionSharedPtr> &EdgeExp,
515  Array<OneD, NekDouble> &edgePhys,
516  const StdRegions::VarCoeffMap &varcoeffs,
517  Array<OneD, NekDouble> &outarray)
518  {
519  bool mmf = (varcoeffs.find(StdRegions::eVarCoeffMF1x) !=
520  varcoeffs.end());
521  int i, j, n;
522  int nquad_e = EdgeExp[edge]->GetNumPoints(0);
523  int order_e = EdgeExp[edge]->GetNcoeffs();
524  int coordim = mmf ? 2 : GetCoordim();
525  int ncoeffs = GetNcoeffs();
526 
527  Array<OneD, NekDouble> inval (nquad_e);
528  Array<OneD, NekDouble> outcoeff(order_e);
529  Array<OneD, NekDouble> tmpcoeff(ncoeffs);
530 
531  const Array<OneD, const Array<OneD, NekDouble> > &normals
532  = GetEdgeNormal(edge);
533 
534  Array<OneD,unsigned int> emap;
535  Array<OneD,int> sign;
536 
537  StdRegions::Orientation edgedir = GetEorient(edge);
538 
539  DNekVec Coeffs (ncoeffs,outarray,eWrapper);
540  DNekVec Tmpcoeff(ncoeffs,tmpcoeff,eWrapper);
541 
542  GetEdgeToElementMap(edge,edgedir,emap,sign);
543 
544  StdRegions::MatrixType DerivType[3] =
545  {
546  StdRegions::eWeakDeriv0,
547  StdRegions::eWeakDeriv1,
548  StdRegions::eWeakDeriv2
549  };
550 
551  StdRegions::VarCoeffType VarCoeff[3] =
552  {
553  StdRegions::eVarCoeffD00,
554  StdRegions::eVarCoeffD11,
555  StdRegions::eVarCoeffD22
556  };
557 
558  StdRegions::VarCoeffMap::const_iterator x;
559  /// @TODO: What direction to use here??
560  if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
561  {
562  Array<OneD, NekDouble> work(nquad_e);
563  GetPhysEdgeVarCoeffsFromElement(
564  edge, EdgeExp[edge], x->second, work);
565  Vmath::Vmul(nquad_e, work, 1, edgePhys, 1, edgePhys, 1);
566  }
567 
568  //================================================================
569  // Add F = \tau <phi_i,in_phys>
570  // Fill edge and take inner product
571  EdgeExp[edge]->IProductWRTBase(edgePhys, outcoeff);
572  // add data to out array
573  for(i = 0; i < order_e; ++i)
574  {
575  outarray[emap[i]] += sign[i] * tau * outcoeff[i];
576  }
577  //================================================================
578 
579  //===============================================================
580  // Add -\sum_i D_i^T M^{-1} G_i + E_i M^{-1} G_i =
581  // \sum_i D_i M^{-1} G_i term
582 
583  // Two independent direction
584  DNekScalMatSharedPtr invMass;
585  for(n = 0; n < coordim; ++n)
586  {
587  if (mmf)
588  {
589  StdRegions::VarCoeffMap Weight;
590  Weight[StdRegions::eVarCoeffMass] = v_GetMFMag(n,varcoeffs);
591 
592  MatrixKey invMasskey( StdRegions::eInvMass,
593  DetShapeType(), *this,
594  StdRegions::NullConstFactorMap,
595  Weight);
596 
597  invMass = GetLocMatrix(invMasskey);
598 
599  Array<OneD, NekDouble> ncdotMF_e =
600  v_GetnEdgecdotMF(n, edge, EdgeExp[edge], normals,
601  varcoeffs);
602 
603  Vmath::Vmul(nquad_e, ncdotMF_e, 1, edgePhys, 1, inval, 1);
604  }
605  else
606  {
607  Vmath::Vmul(nquad_e, normals[n], 1, edgePhys, 1, inval, 1);
608  invMass = GetLocMatrix(StdRegions::eInvMass);
609  }
610 
611  if (m_negatedNormals[edge])
612  {
613  Vmath::Neg(nquad_e, inval, 1);
614  }
615 
616  // Multiply by variable coefficient
617  /// @TODO: Document this (probably not needed)
618 // StdRegions::VarCoeffMap::const_iterator x;
619 // if ((x = varcoeffs.find(VarCoeff[n])) != varcoeffs.end())
620 // {
621 // GetPhysEdgeVarCoeffsFromElement(edge,EdgeExp[edge],x->second,varcoeff_work);
622 // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[edge]->GetPhys(),1,EdgeExp[edge]->UpdatePhys(),1);
623 // }
624 
625  EdgeExp[edge]->IProductWRTBase(inval, outcoeff);
626 
627  // M^{-1} G
628  for(i = 0; i < ncoeffs; ++i)
629  {
630  tmpcoeff[i] = 0;
631  for(j = 0; j < order_e; ++j)
632  {
633  tmpcoeff[i] += (*invMass)(i,emap[j]) * sign[j]
634  * outcoeff[j];
635  }
636  }
637 
638  if (mmf)
639  {
640  StdRegions::VarCoeffMap VarCoeffDirDeriv;
641  VarCoeffDirDeriv[StdRegions::eVarCoeffMF] =
642  v_GetMF(n,coordim,varcoeffs);
643  VarCoeffDirDeriv[StdRegions::eVarCoeffMFDiv] =
644  v_GetMFDiv(n,varcoeffs);
645 
646  MatrixKey Dmatkey( StdRegions::eWeakDirectionalDeriv,
647  DetShapeType(), *this,
648  StdRegions::NullConstFactorMap,
649  VarCoeffDirDeriv);
650 
651  DNekScalMat &Dmat = *GetLocMatrix(Dmatkey);
652 
653  Coeffs = Coeffs + Dmat*Tmpcoeff;
654  }
655  else
656  {
657  if(varcoeffs.find(VarCoeff[n]) != varcoeffs.end())
658  {
659  MatrixKey mkey(DerivType[n], DetShapeType(), *this,
660  StdRegions::NullConstFactorMap, varcoeffs);
661 
662  DNekScalMat &Dmat = *GetLocMatrix(mkey);
663  Coeffs = Coeffs + Dmat*Tmpcoeff;
664  }
665  else
666  {
667  DNekScalMat &Dmat = *GetLocMatrix(DerivType[n]);
668  Coeffs = Coeffs + Dmat*Tmpcoeff;
669  }
670  }
671  }
672  }
673 
674 
675  /**
676  * Extracts the variable coefficients along an edge
677  */
678  void Expansion2D::GetPhysEdgeVarCoeffsFromElement(
679  const int edge,
680  ExpansionSharedPtr &EdgeExp,
681  const Array<OneD, const NekDouble> &varcoeff,
682  Array<OneD,NekDouble> &outarray)
683  {
684  Array<OneD, NekDouble> tmp(GetNcoeffs());
685  Array<OneD, NekDouble> edgetmp(EdgeExp->GetNcoeffs());
686 
687  // FwdTrans varcoeffs
688  FwdTrans(varcoeff, tmp);
689 
690  // Map to edge
691  Array<OneD,unsigned int> emap;
692  Array<OneD, int> sign;
693  StdRegions::Orientation edgedir = GetEorient(edge);
694  GetEdgeToElementMap(edge,edgedir,emap,sign);
695 
696  for (unsigned int i = 0; i < EdgeExp->GetNcoeffs(); ++i)
697  {
698  edgetmp[i] = tmp[emap[i]];
699  }
700 
701  // BwdTrans
702  EdgeExp->BwdTrans(edgetmp, outarray);
703  }
704 
705 
706  /**
707  * Computes matrices needed for the HDG formulation. References to
708  * equations relate to the following paper:
709  * R. M. Kirby, S. J. Sherwin, B. Cockburn, To CG or to HDG: A
710  * Comparative Study, J. Sci. Comp P1-30
711  * DOI 10.1007/s10915-011-9501-7
712  */
713  DNekMatSharedPtr Expansion2D::v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
714  {
715  DNekMatSharedPtr returnval;
716 
717  switch(mkey.GetMatrixType())
718  {
719  // (Z^e)^{-1} (Eqn. 33, P22)
720  case StdRegions::eHybridDGHelmholtz:
721  {
722  ASSERTL1(IsBoundaryInteriorExpansion(),
723  "HybridDGHelmholtz matrix not set up "
724  "for non boundary-interior expansions");
725 
726  int i,j,k;
727  NekDouble lambdaval = mkey.GetConstFactor(StdRegions::eFactorLambda);
728  NekDouble tau = mkey.GetConstFactor(StdRegions::eFactorTau);
729  int ncoeffs = GetNcoeffs();
730  int nedges = GetNedges();
731  int shapedim = 2;
732  const StdRegions::VarCoeffMap &varcoeffs
733  = mkey.GetVarCoeffs();
734  bool mmf =
735  (varcoeffs.find(StdRegions::eVarCoeffMF1x) !=
736  varcoeffs.end());
737 
738  Array<OneD,unsigned int> emap;
739  Array<OneD,int> sign;
740  StdRegions::Orientation edgedir = StdRegions::eForwards;
741  ExpansionSharedPtr EdgeExp;
742  ExpansionSharedPtr EdgeExp2;
743 
744  int order_e, coordim = GetCoordim();
745  DNekScalMat &invMass = *GetLocMatrix(StdRegions::eInvMass);
746  StdRegions::MatrixType DerivType[3] = {StdRegions::eWeakDeriv0,
747  StdRegions::eWeakDeriv1,
748  StdRegions::eWeakDeriv2};
749  DNekMat LocMat(ncoeffs,ncoeffs);
750 
751  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(ncoeffs,ncoeffs);
752  DNekMat &Mat = *returnval;
753  Vmath::Zero(ncoeffs*ncoeffs,Mat.GetPtr(),1);
754 
755  StdRegions::VarCoeffType Coeffs[3] = {StdRegions::eVarCoeffD00,
756  StdRegions::eVarCoeffD11,
757  StdRegions::eVarCoeffD22};
758 
759  StdRegions::VarCoeffMap::const_iterator x;
760 
761  for(i=0; i < coordim; ++i)
762  {
763  if (mmf)
764  {
765  if(i < shapedim)
766  {
767  StdRegions::VarCoeffMap VarCoeffDirDeriv;
768  VarCoeffDirDeriv[StdRegions::eVarCoeffMF] =
769  v_GetMF(i,shapedim,varcoeffs);
770  VarCoeffDirDeriv[StdRegions::eVarCoeffMFDiv] =
771  v_GetMFDiv(i,varcoeffs);
772 
773  MatrixKey Dmatkey(StdRegions::
774  eWeakDirectionalDeriv,
775  DetShapeType(), *this,
776  StdRegions::
777  NullConstFactorMap,
778  VarCoeffDirDeriv);
779 
780  DNekScalMat &Dmat = *GetLocMatrix(Dmatkey);
781 
782  StdRegions::VarCoeffMap Weight;
783  Weight[StdRegions::eVarCoeffMass] =
784  v_GetMFMag(i,mkey.GetVarCoeffs());
785 
786  MatrixKey invMasskey( StdRegions::eInvMass,
787  DetShapeType(), *this,
788  StdRegions::
789  NullConstFactorMap,
790  Weight);
791 
792  DNekScalMat &invMass =
793  *GetLocMatrix(invMasskey);
794 
795  Mat = Mat + Dmat*invMass*Transpose(Dmat);
796  }
797  }
798  else if(mkey.HasVarCoeff(Coeffs[i]))
799  {
800  MatrixKey DmatkeyL( DerivType[i],
801  DetShapeType(), *this,
802  StdRegions::NullConstFactorMap,
803  mkey.GetVarCoeffAsMap(
804  Coeffs[i]));
805 
806  MatrixKey DmatkeyR( DerivType[i],
807  DetShapeType(), *this);
808 
809  DNekScalMat &DmatL = *GetLocMatrix(DmatkeyL);
810  DNekScalMat &DmatR = *GetLocMatrix(DmatkeyR);
811  Mat = Mat + DmatL*invMass*Transpose(DmatR);
812  }
813  else
814  {
815  DNekScalMat &Dmat = *GetLocMatrix(DerivType[i]);
816  Mat = Mat + Dmat*invMass*Transpose(Dmat);
817  }
818 
819  }
820 
821  // Add Mass Matrix Contribution for Helmholtz problem
822  DNekScalMat &Mass = *GetLocMatrix(StdRegions::eMass);
823  Mat = Mat + lambdaval*Mass;
824 
825  // Add tau*E_l using elemental mass matrices on each edge
826  for(i = 0; i < nedges; ++i)
827  {
828  EdgeExp = GetEdgeExp(i);
829  EdgeExp2 = GetEdgeExp(i);
830  order_e = EdgeExp->GetNcoeffs();
831  int nq = EdgeExp->GetNumPoints(0);
832  GetEdgeToElementMap(i,edgedir,emap,sign);
833 
834  // @TODO: Document
835  StdRegions::VarCoeffMap edgeVarCoeffs;
836  if (mkey.HasVarCoeff(StdRegions::eVarCoeffD00))
837  {
838  Array<OneD, NekDouble> mu(nq);
839  GetPhysEdgeVarCoeffsFromElement(i, EdgeExp2, mkey.GetVarCoeff(StdRegions::eVarCoeffD00), mu);
840  edgeVarCoeffs[StdRegions::eVarCoeffMass] = mu;
841  }
842  DNekScalMat &eMass = *EdgeExp->GetLocMatrix(StdRegions::eMass, StdRegions::NullConstFactorMap, edgeVarCoeffs);
843  //DNekScalMat &eMass = *EdgeExp->GetLocMatrix(StdRegions::eMass);
844 
845  for(j = 0; j < order_e; ++j)
846  {
847  for(k = 0; k < order_e; ++k)
848  {
849  Mat(emap[j],emap[k]) = Mat(emap[j],emap[k]) + tau*sign[j]*sign[k]*eMass(j,k);
850  }
851  }
852  }
853  }
854  break;
855  // U^e (P22)
856  case StdRegions::eHybridDGLamToU:
857  {
858  int i,j,k;
859  int nbndry = NumDGBndryCoeffs();
860  int ncoeffs = GetNcoeffs();
861  int nedges = GetNedges();
862  NekDouble tau = mkey.GetConstFactor(StdRegions::eFactorTau);
863 
864  Array<OneD,NekDouble> lambda(nbndry);
865  DNekVec Lambda(nbndry,lambda,eWrapper);
866  Array<OneD,NekDouble> ulam(ncoeffs);
867  DNekVec Ulam(ncoeffs,ulam,eWrapper);
868  Array<OneD,NekDouble> f(ncoeffs);
869  DNekVec F(ncoeffs,f,eWrapper);
870 
871  Array<OneD, ExpansionSharedPtr> EdgeExp(nedges);
872  // declare matrix space
873  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(ncoeffs,nbndry);
874  DNekMat &Umat = *returnval;
875 
876  // Z^e matrix
877  MatrixKey newkey(StdRegions::eInvHybridDGHelmholtz, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
878  DNekScalMat &invHmat = *GetLocMatrix(newkey);
879 
880  Array<OneD,unsigned int> emap;
881  Array<OneD,int> sign;
882 
883  for(i = 0; i < nedges; ++i)
884  {
885  EdgeExp[i] = GetEdgeExp(i);
886  }
887 
888  // for each degree of freedom of the lambda space
889  // calculate Umat entry
890  // Generate Lambda to U_lambda matrix
891  for(j = 0; j < nbndry; ++j)
892  {
893  // standard basis vectors e_j
894  Vmath::Zero(nbndry,&lambda[0],1);
895  Vmath::Zero(ncoeffs,&f[0],1);
896  lambda[j] = 1.0;
897 
898  SetTraceToGeomOrientation(EdgeExp,lambda);
899 
900  // Compute F = [I D_1 M^{-1} D_2 M^{-1}] C e_j
901  AddHDGHelmholtzTraceTerms(
902  tau, lambda, EdgeExp, mkey.GetVarCoeffs(), f);
903 
904  // Compute U^e_j
905  Ulam = invHmat*F; // generate Ulam from lambda
906 
907  // fill column of matrix
908  for(k = 0; k < ncoeffs; ++k)
909  {
910  Umat(k,j) = Ulam[k];
911  }
912  }
913  }
914  break;
915  // Q_0, Q_1, Q_2 matrices (P23)
916  // Each are a product of a row of Eqn 32 with the C matrix.
917  // Rather than explicitly computing all of Eqn 32, we note each
918  // row is almost a multiple of U^e, so use that as our starting
919  // point.
920  case StdRegions::eHybridDGLamToQ0:
921  case StdRegions::eHybridDGLamToQ1:
922  case StdRegions::eHybridDGLamToQ2:
923  {
924  int i = 0;
925  int j = 0;
926  int k = 0;
927  int dir = 0;
928  int nbndry = NumDGBndryCoeffs();
929  int ncoeffs = GetNcoeffs();
930  int nedges = GetNedges();
931  int shapedim = 2;
932 
933  Array<OneD,NekDouble> lambda(nbndry);
934  DNekVec Lambda(nbndry,lambda,eWrapper);
935  Array<OneD, ExpansionSharedPtr> EdgeExp(nedges);
936 
937  Array<OneD,NekDouble> ulam(ncoeffs);
938  DNekVec Ulam(ncoeffs,ulam,eWrapper);
939  Array<OneD,NekDouble> f(ncoeffs);
940  DNekVec F(ncoeffs,f,eWrapper);
941 
942  // declare matrix space
943  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(ncoeffs,nbndry);
944  DNekMat &Qmat = *returnval;
945 
946  // Lambda to U matrix
947  MatrixKey lamToUkey(StdRegions::eHybridDGLamToU, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
948  DNekScalMat &lamToU = *GetLocMatrix(lamToUkey);
949 
950  // Inverse mass matrix
951  DNekScalMat &invMass = *GetLocMatrix(StdRegions::eInvMass);
952 
953 
954  for(i = 0; i < nedges; ++i)
955  {
956  EdgeExp[i] = GetEdgeExp(i);
957  }
958 
959  //Weak Derivative matrix
960  DNekScalMatSharedPtr Dmat;
961  switch(mkey.GetMatrixType())
962  {
963  case StdRegions::eHybridDGLamToQ0:
964  dir = 0;
965  Dmat = GetLocMatrix(StdRegions::eWeakDeriv0);
966  break;
967  case StdRegions::eHybridDGLamToQ1:
968  dir = 1;
969  Dmat = GetLocMatrix(StdRegions::eWeakDeriv1);
970  break;
971  case StdRegions::eHybridDGLamToQ2:
972  dir = 2;
973  Dmat = GetLocMatrix(StdRegions::eWeakDeriv2);
974  break;
975  default:
976  ASSERTL0(false,"Direction not known");
977  break;
978  }
979 
980  const StdRegions::VarCoeffMap &varcoeffs =
981  mkey.GetVarCoeffs();
982  if (varcoeffs.find(StdRegions::eVarCoeffMF1x) !=
983  varcoeffs.end())
984  {
985  StdRegions::VarCoeffMap VarCoeffDirDeriv;
986  VarCoeffDirDeriv[StdRegions::eVarCoeffMF] =
987  v_GetMF(dir,shapedim,varcoeffs);
988  VarCoeffDirDeriv[StdRegions::eVarCoeffMFDiv] =
989  v_GetMFDiv(dir,varcoeffs);
990 
991  MatrixKey Dmatkey( StdRegions::eWeakDirectionalDeriv,
992  DetShapeType(), *this,
993  StdRegions::NullConstFactorMap,
994  VarCoeffDirDeriv);
995 
996  Dmat = GetLocMatrix(Dmatkey);
997 
998  StdRegions::VarCoeffMap Weight;
999  Weight[StdRegions::eVarCoeffMass] =
1000  v_GetMFMag(dir,mkey.GetVarCoeffs());
1001 
1002  MatrixKey invMasskey( StdRegions::eInvMass,
1003  DetShapeType(), *this,
1004  StdRegions::NullConstFactorMap,
1005  Weight);
1006 
1007  invMass = *GetLocMatrix(invMasskey);
1008  }
1009  else
1010  {
1011  StdRegions::MatrixType DerivType[3] = {
1012  StdRegions::eWeakDeriv0,
1013  StdRegions::eWeakDeriv1,
1014  StdRegions::eWeakDeriv2};
1015 
1016  Dmat = GetLocMatrix(DerivType[dir]);
1017 
1018  MatrixKey invMasskey( StdRegions::eInvMass,
1019  DetShapeType(), *this);
1020  invMass = *GetLocMatrix(invMasskey);
1021  }
1022 
1023  // for each degree of freedom of the lambda space
1024  // calculate Qmat entry
1025  // Generate Lambda to Q_lambda matrix
1026  for(j = 0; j < nbndry; ++j)
1027  {
1028  Vmath::Zero(nbndry,&lambda[0],1);
1029  lambda[j] = 1.0;
1030 
1031  // for lambda[j] = 1 this is the solution to ulam
1032  for(k = 0; k < ncoeffs; ++k)
1033  {
1034  Ulam[k] = lamToU(k,j);
1035  }
1036 
1037  // -D^T ulam
1038  Vmath::Neg(ncoeffs,&ulam[0],1);
1039  F = Transpose(*Dmat)*Ulam;
1040 
1041  SetTraceToGeomOrientation(EdgeExp,lambda);
1042 
1043  // Add the C terms resulting from the I's on the
1044  // diagonals of Eqn 32
1045  AddNormTraceInt(dir,lambda,EdgeExp,f,mkey.GetVarCoeffs());
1046 
1047  // finally multiply by inverse mass matrix
1048  Ulam = invMass*F;
1049 
1050  // fill column of matrix (Qmat is in column major format)
1051  Vmath::Vcopy(ncoeffs,&ulam[0],1,&(Qmat.GetPtr())[0]+j*ncoeffs,1);
1052  }
1053  }
1054  break;
1055  // Matrix K (P23)
1056  case StdRegions::eHybridDGHelmBndLam:
1057  {
1058  int i,j,e,cnt;
1059  int order_e, nquad_e;
1060  int nbndry = NumDGBndryCoeffs();
1061  int coordim = GetCoordim();
1062  int nedges = GetNedges();
1063  NekDouble tau = mkey.GetConstFactor(StdRegions::eFactorTau);
1064  StdRegions::VarCoeffMap::const_iterator x;
1065  const StdRegions::VarCoeffMap &varcoeffs
1066  = mkey.GetVarCoeffs();
1067  bool mmf =
1068  (varcoeffs.find(StdRegions::eVarCoeffMF1x) !=
1069  varcoeffs.end());
1070 
1071  Array<OneD,NekDouble> work, varcoeff_work;
1072  Array<OneD,const Array<OneD, NekDouble> > normals;
1073  Array<OneD, ExpansionSharedPtr> EdgeExp(nedges);
1074  Array<OneD, NekDouble> lam(nbndry);
1075 
1076  Array<OneD,unsigned int> emap;
1077  Array<OneD, int> sign;
1078  StdRegions::Orientation edgedir;
1079 
1080  // declare matrix space
1081  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(nbndry, nbndry);
1082  DNekMat &BndMat = *returnval;
1083 
1084  DNekScalMatSharedPtr LamToQ[3];
1085 
1086  // Matrix to map Lambda to U
1087  MatrixKey LamToUkey(StdRegions::eHybridDGLamToU, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1088  DNekScalMat &LamToU = *GetLocMatrix(LamToUkey);
1089 
1090  // Matrix to map Lambda to Q0
1091  MatrixKey LamToQ0key(StdRegions::eHybridDGLamToQ0, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1092  LamToQ[0] = GetLocMatrix(LamToQ0key);
1093 
1094  // Matrix to map Lambda to Q1
1095  MatrixKey LamToQ1key(StdRegions::eHybridDGLamToQ1, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1096  LamToQ[1] = GetLocMatrix(LamToQ1key);
1097 
1098  // Matrix to map Lambda to Q2 for 3D coordinates
1099  if (coordim == 3)
1100  {
1101  MatrixKey LamToQ2key(StdRegions::eHybridDGLamToQ2, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1102  LamToQ[2] = GetLocMatrix(LamToQ2key);
1103  }
1104 
1105  // Set up edge segment expansions from local geom info
1106  for(i = 0; i < nedges; ++i)
1107  {
1108  EdgeExp[i] = GetEdgeExp(i);
1109  }
1110 
1111  // Set up matrix derived from <mu, Q_lam.n - \tau (U_lam - Lam) >
1112  for(i = 0; i < nbndry; ++i)
1113  {
1114  cnt = 0;
1115 
1116  Vmath::Zero(nbndry,lam,1);
1117  lam[i] = 1.0;
1118  SetTraceToGeomOrientation(EdgeExp,lam);
1119 
1120  for(e = 0; e < nedges; ++e)
1121  {
1122  order_e = EdgeExp[e]->GetNcoeffs();
1123  nquad_e = EdgeExp[e]->GetNumPoints(0);
1124 
1125  normals = GetEdgeNormal(e);
1126  edgedir = GetEorient(e);
1127 
1128  work = Array<OneD,NekDouble>(nquad_e);
1129  varcoeff_work = Array<OneD, NekDouble>(nquad_e);
1130 
1131  GetEdgeToElementMap(e,edgedir,emap,sign);
1132 
1133 
1134  StdRegions::VarCoeffType VarCoeff[3] = {StdRegions::eVarCoeffD00,
1135  StdRegions::eVarCoeffD11,
1136  StdRegions::eVarCoeffD22};
1137 
1138  // Q0 * n0 (BQ_0 terms)
1139  Array<OneD, NekDouble> edgeCoeffs(order_e);
1140  Array<OneD, NekDouble> edgePhys (nquad_e);
1141  for(j = 0; j < order_e; ++j)
1142  {
1143  edgeCoeffs[j] = sign[j]*(*LamToQ[0])(emap[j],i);
1144  }
1145 
1146  EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
1147 // @TODO Var coeffs
1148  // Multiply by variable coefficient
1149 // if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
1150 // {
1151 // GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
1152 // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
1153 // }
1154  if (mmf)
1155  {
1156  Array<OneD, NekDouble> ncdotMF =
1157  v_GetnEdgecdotMF(0, e, EdgeExp[e],
1158  normals, varcoeffs);
1159  Vmath::Vmul(nquad_e, ncdotMF, 1,
1160  edgePhys, 1,
1161  work, 1);
1162  }
1163  else
1164  {
1165  Vmath::Vmul(nquad_e, normals[0], 1,
1166  edgePhys, 1,
1167  work, 1);
1168  }
1169 
1170  // Q1 * n1 (BQ_1 terms)
1171  for(j = 0; j < order_e; ++j)
1172  {
1173  edgeCoeffs[j] = sign[j]*(*LamToQ[1])(emap[j],i);
1174  }
1175 
1176  EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
1177 
1178 // @TODO var coeffs
1179  // Multiply by variable coefficients
1180 // if ((x = varcoeffs.find(VarCoeff[1])) != varcoeffs.end())
1181 // {
1182 // GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
1183 // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
1184 // }
1185 
1186  if (mmf)
1187  {
1188  Array<OneD, NekDouble> ncdotMF =
1189  v_GetnEdgecdotMF(1, e, EdgeExp[e],
1190  normals, varcoeffs);
1191  Vmath::Vvtvp(nquad_e, ncdotMF, 1,
1192  edgePhys, 1,
1193  work, 1,
1194  work, 1);
1195  }
1196  else
1197  {
1198  Vmath::Vvtvp(nquad_e, normals[1], 1,
1199  edgePhys, 1,
1200  work, 1,
1201  work, 1);
1202  }
1203 
1204  // Q2 * n2 (BQ_2 terms)
1205  if (coordim == 3)
1206  {
1207  for(j = 0; j < order_e; ++j)
1208  {
1209  edgeCoeffs[j] = sign[j]*(*LamToQ[2])(emap[j],i);
1210  }
1211 
1212  EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
1213 // @TODO var coeffs
1214  // Multiply by variable coefficients
1215 // if ((x = varcoeffs.find(VarCoeff[2])) != varcoeffs.end())
1216 // {
1217 // GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
1218 // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
1219 // }
1220 
1221  Vmath::Vvtvp(nquad_e, normals[2], 1, edgePhys, 1,
1222  work, 1, work, 1);
1223  }
1224 
1225  if (m_negatedNormals[e])
1226  {
1227  Vmath::Neg(nquad_e, work, 1);
1228  }
1229 
1230  // - tau (ulam - lam)
1231  // Corresponds to the G and BU terms.
1232  for(j = 0; j < order_e; ++j)
1233  {
1234  edgeCoeffs[j] = sign[j]*LamToU(emap[j],i) - lam[cnt+j];
1235  }
1236 
1237  EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
1238 
1239  // Multiply by variable coefficients
1240  if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
1241  {
1242  GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
1243  Vmath::Vmul(nquad_e,varcoeff_work,1,edgePhys,1,edgePhys,1);
1244  }
1245 
1246  Vmath::Svtvp(nquad_e,-tau,edgePhys,1,
1247  work,1,work,1);
1248 /// TODO: Add variable coeffs
1249  EdgeExp[e]->IProductWRTBase(work, edgeCoeffs);
1250 
1251  EdgeExp[e]->SetCoeffsToOrientation(edgeCoeffs, edgedir);
1252 
1253  for(j = 0; j < order_e; ++j)
1254  {
1255  BndMat(cnt+j,i) = edgeCoeffs[j];
1256  }
1257 
1258  cnt += order_e;
1259  }
1260  }
1261  }
1262  break;
1263  //HDG postprocessing
1264  case StdRegions::eInvLaplacianWithUnityMean:
1265  {
1266  MatrixKey lapkey(StdRegions::eLaplacian, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1267  DNekScalMat &LapMat = *GetLocMatrix(lapkey);
1268 
1269  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(LapMat.GetRows(),LapMat.GetColumns());
1270  DNekMatSharedPtr lmat = returnval;
1271 
1272  (*lmat) = LapMat;
1273 
1274  // replace first column with inner product wrt 1
1275  int nq = GetTotPoints();
1276  Array<OneD, NekDouble> tmp(nq);
1277  Array<OneD, NekDouble> outarray(m_ncoeffs);
1278  Vmath::Fill(nq,1.0,tmp,1);
1279  IProductWRTBase(tmp, outarray);
1280 
1281  Vmath::Vcopy(m_ncoeffs,&outarray[0],1,
1282  &(lmat->GetPtr())[0],1);
1283 
1284  lmat->Invert();
1285  }
1286  break;
1287  default:
1288  ASSERTL0(false,"This matrix type cannot be generated from this class");
1289  break;
1290  }
1291 
1292  return returnval;
1293  }
1294 
1295  //Evaluate Coefficients of weak deriviative in the direction dir
1296  //given the input coefficicents incoeffs and the imposed
1297  //boundary values in EdgeExp (which will have its phys space updated);
1298  void Expansion2D::v_DGDeriv(
1299  int dir,
1300  const Array<OneD, const NekDouble> &incoeffs,
1301  Array<OneD, ExpansionSharedPtr> &EdgeExp,
1302  Array<OneD, Array<OneD, NekDouble> > &edgeCoeffs,
1303  Array<OneD, NekDouble> &out_d)
1304  {
1305  StdRegions::MatrixType DerivType[3] = {StdRegions::eWeakDeriv0,
1306  StdRegions::eWeakDeriv1,
1307  StdRegions::eWeakDeriv2};
1308 
1309  int ncoeffs = GetNcoeffs();
1310 
1311  DNekScalMat &InvMass = *GetLocMatrix(StdRegions::eInvMass);
1312  DNekScalMat &Dmat = *GetLocMatrix(DerivType[dir]);
1313 
1314  Array<OneD, NekDouble> coeffs = incoeffs;
1315  DNekVec Coeffs (ncoeffs,coeffs, eWrapper);
1316 
1317  Coeffs = Transpose(Dmat)*Coeffs;
1318  Vmath::Neg(ncoeffs, coeffs,1);
1319 
1320  // Add the boundary integral including the relevant part of
1321  // the normal
1322  AddNormTraceInt(dir, EdgeExp, edgeCoeffs, coeffs);
1323 
1324  DNekVec Out_d (ncoeffs,out_d,eWrapper);
1325 
1326  Out_d = InvMass*Coeffs;
1327  }
1328 
1329  enum BndToLocMatrixMapType
1330  {
1331  eBndToFullMatrixCG,
1332  eBndToBndMatrixCG,
1333  eBndToTraceMatrixDG
1334  };
1335 
1336  void Expansion2D::v_AddRobinMassMatrix(const int edge, const Array<OneD, const NekDouble > &primCoeffs, DNekMatSharedPtr &inoutmat)
1337  {
1338  ASSERTL1(IsBoundaryInteriorExpansion(),
1339  "Not set up for non boundary-interior expansions");
1340  ASSERTL1(inoutmat->GetRows() == inoutmat->GetColumns(),
1341  "Assuming that input matrix was square");
1342  int i,j;
1343  int id1,id2;
1344  ExpansionSharedPtr edgeExp = m_edgeExp[edge].lock();
1345  int order_e = edgeExp->GetNcoeffs();
1346 
1347  Array<OneD,unsigned int> map;
1348  Array<OneD,int> sign;
1349 
1350  StdRegions::VarCoeffMap varcoeffs;
1351  varcoeffs[StdRegions::eVarCoeffMass] = primCoeffs;
1352 
1353  LocalRegions::MatrixKey mkey(StdRegions::eMass,LibUtilities::eSegment, *edgeExp, StdRegions::NullConstFactorMap, varcoeffs);
1354  DNekScalMat &edgemat = *edgeExp->GetLocMatrix(mkey);
1355 
1356  // Now need to identify a map which takes the local edge
1357  // mass matrix to the matrix stored in inoutmat;
1358  // This can currently be deduced from the size of the matrix
1359 
1360  // - if inoutmat.m_rows() == v_NCoeffs() it is a full
1361  // matrix system
1362 
1363  // - if inoutmat.m_rows() == v_NumBndCoeffs() it is a
1364  // boundary CG system
1365 
1366  // - if inoutmat.m_rows() == v_NumDGBndCoeffs() it is a
1367  // trace DG system
1368  int rows = inoutmat->GetRows();
1369 
1370  if (rows == GetNcoeffs())
1371  {
1372  GetEdgeToElementMap(edge,v_GetEorient(edge),map,sign);
1373  }
1374  else if(rows == NumBndryCoeffs())
1375  {
1376  int nbndry = NumBndryCoeffs();
1377  Array<OneD,unsigned int> bmap(nbndry);
1378 
1379  GetEdgeToElementMap(edge,v_GetEorient(edge),map,sign);
1380 
1381  GetBoundaryMap(bmap);
1382 
1383  for(i = 0; i < order_e; ++i)
1384  {
1385  for(j = 0; j < nbndry; ++j)
1386  {
1387  if(map[i] == bmap[j])
1388  {
1389  map[i] = j;
1390  break;
1391  }
1392  }
1393  ASSERTL1(j != nbndry,"Did not find number in map");
1394  }
1395  }
1396  else if (rows == NumDGBndryCoeffs())
1397  {
1398  // possibly this should be a separate method
1399  int cnt = 0;
1400  map = Array<OneD, unsigned int> (order_e);
1401  sign = Array<OneD, int> (order_e,1);
1402 
1403  for(i = 0; i < edge; ++i)
1404  {
1405  cnt += GetEdgeNcoeffs(i);
1406  }
1407 
1408  for(i = 0; i < order_e; ++i)
1409  {
1410  map[i] = cnt++;
1411  }
1412  // check for mapping reversal
1413  if(GetEorient(edge) == StdRegions::eBackwards)
1414  {
1415  switch(edgeExp->GetBasis(0)->GetBasisType())
1416  {
1417  case LibUtilities::eGauss_Lagrange:
1418  reverse( map.get() , map.get()+order_e);
1419  break;
1420  case LibUtilities::eGLL_Lagrange:
1421  reverse( map.get() , map.get()+order_e);
1422  break;
1423  case LibUtilities::eModified_A:
1424  {
1425  swap(map[0],map[1]);
1426  for(i = 3; i < order_e; i+=2)
1427  {
1428  sign[i] = -1;
1429  }
1430  }
1431  break;
1432  default:
1433  ASSERTL0(false,"Edge boundary type not valid for this method");
1434  }
1435  }
1436  }
1437  else
1438  {
1439  ASSERTL0(false,"Could not identify matrix type from dimension");
1440  }
1441 
1442  for(i = 0; i < order_e; ++i)
1443  {
1444  id1 = map[i];
1445  for(j = 0; j < order_e; ++j)
1446  {
1447  id2 = map[j];
1448  (*inoutmat)(id1,id2) += edgemat(i,j)*sign[i]*sign[j];
1449  }
1450  }
1451  }
1452 
1453  /**
1454  * Given an edge and vector of element coefficients:
1455  * - maps those elemental coefficients corresponding to the edge into
1456  * an edge-vector.
1457  * - resets the element coefficients
1458  * - multiplies the edge vector by the edge mass matrix
1459  * - maps the edge coefficients back onto the elemental coefficients
1460  */
1461  void Expansion2D::v_AddRobinEdgeContribution(const int edgeid, const Array<OneD, const NekDouble> &primCoeffs, Array<OneD, NekDouble> &coeffs)
1462  {
1463  ASSERTL1(IsBoundaryInteriorExpansion(),
1464  "Not set up for non boundary-interior expansions");
1465  int i;
1466  ExpansionSharedPtr edgeExp = m_edgeExp[edgeid].lock();
1467  int order_e = edgeExp->GetNcoeffs();
1468 
1469  Array<OneD,unsigned int> map;
1470  Array<OneD,int> sign;
1471 
1472  StdRegions::VarCoeffMap varcoeffs;
1473  varcoeffs[StdRegions::eVarCoeffMass] = primCoeffs;
1474 
1475  LocalRegions::MatrixKey mkey(StdRegions::eMass,LibUtilities::eSegment, *edgeExp, StdRegions::NullConstFactorMap, varcoeffs);
1476  DNekScalMat &edgemat = *edgeExp->GetLocMatrix(mkey);
1477 
1478  NekVector<NekDouble> vEdgeCoeffs (order_e);
1479 
1480  GetEdgeToElementMap(edgeid,v_GetEorient(edgeid),map,sign);
1481 
1482  for (i = 0; i < order_e; ++i)
1483  {
1484  vEdgeCoeffs[i] = coeffs[map[i]]*sign[i];
1485  }
1486  Vmath::Zero(GetNcoeffs(), coeffs, 1);
1487 
1488  vEdgeCoeffs = edgemat * vEdgeCoeffs;
1489 
1490  for (i = 0; i < order_e; ++i)
1491  {
1492  coeffs[map[i]] = vEdgeCoeffs[i]*sign[i];
1493  }
1494  }
1495 
1496  DNekMatSharedPtr Expansion2D::v_BuildVertexMatrix(
1497  const DNekScalMatSharedPtr &r_bnd)
1498  {
1499  MatrixStorage storage = eFULL;
1500  DNekMatSharedPtr m_vertexmatrix;
1501 
1502  int nVerts, vid1, vid2, vMap1, vMap2;
1503  NekDouble VertexValue;
1504 
1505  nVerts = GetNverts();
1506 
1507  m_vertexmatrix =
1508  MemoryManager<DNekMat>::AllocateSharedPtr(
1509  nVerts, nVerts, 0.0, storage);
1510  DNekMat &VertexMat = (*m_vertexmatrix);
1511 
1512  for (vid1 = 0; vid1 < nVerts; ++vid1)
1513  {
1514  vMap1 = GetVertexMap(vid1);
1515 
1516  for (vid2 = 0; vid2 < nVerts; ++vid2)
1517  {
1518  vMap2 = GetVertexMap(vid2);
1519  VertexValue = (*r_bnd)(vMap1, vMap2);
1520  VertexMat.SetValue(vid1, vid2, VertexValue);
1521  }
1522  }
1523 
1524  return m_vertexmatrix;
1525  }
1526 
1527  Array<OneD, unsigned int> Expansion2D::v_GetEdgeInverseBoundaryMap(
1528  int eid)
1529  {
1530  int n, j;
1531  int nEdgeCoeffs;
1532  int nBndCoeffs = NumBndryCoeffs();
1533 
1534  Array<OneD, unsigned int> bmap(nBndCoeffs);
1535  GetBoundaryMap(bmap);
1536 
1537  // Map from full system to statically condensed system (i.e reverse
1538  // GetBoundaryMap)
1539  map<int, int> invmap;
1540  for (j = 0; j < nBndCoeffs; ++j)
1541  {
1542  invmap[bmap[j]] = j;
1543  }
1544 
1545  // Number of interior edge coefficients
1546  nEdgeCoeffs = GetEdgeNcoeffs(eid) - 2;
1547 
1548  const SpatialDomains::Geometry2DSharedPtr &geom = GetGeom2D();
1549 
1550  Array<OneD, unsigned int> edgemaparray(nEdgeCoeffs);
1551  Array<OneD, unsigned int> maparray (nEdgeCoeffs);
1552  Array<OneD, int> signarray (nEdgeCoeffs, 1);
1553  StdRegions::Orientation eOrient = geom->GetEorient(eid);
1554 
1555  // maparray is the location of the edge within the matrix
1556  GetEdgeInteriorMap(eid, eOrient, maparray, signarray);
1557 
1558  for (n = 0; n < nEdgeCoeffs; ++n)
1559  {
1560  edgemaparray[n] = invmap[maparray[n]];
1561  }
1562 
1563  return edgemaparray;
1564  }
1565 
1566  void Expansion2D::v_SetUpPhysNormals(const int edge)
1567  {
1568  ComputeEdgeNormal(edge);
1569  }
1570 
1571  const StdRegions::NormalVector &Expansion2D::v_GetEdgeNormal(const int edge) const
1572  {
1573  std::map<int, StdRegions::NormalVector>::const_iterator x;
1574  x = m_edgeNormals.find(edge);
1575  ASSERTL0 (x != m_edgeNormals.end(),
1576  "Edge normal not computed.");
1577  return x->second;
1578  }
1579 
1580  const StdRegions::NormalVector &Expansion2D::v_GetSurfaceNormal(
1581  const int id) const
1582  {
1583  return v_GetEdgeNormal(id);
1584  }
1585 
1586  void Expansion2D::v_NegateEdgeNormal(const int edge)
1587  {
1588  m_negatedNormals[edge] = true;
1589  for (int i = 0; i < GetCoordim(); ++i)
1590  {
1591  Vmath::Neg(m_edgeNormals[edge][i].num_elements(),
1592  m_edgeNormals[edge][i], 1);
1593  }
1594  }
1595 
1596  bool Expansion2D::v_EdgeNormalNegated(const int edge)
1597  {
1598  return m_negatedNormals[edge];
1599  }
1600 
1601  void Expansion2D::ReOrientEdgePhysMap(
1602  const int nvert,
1603  const StdRegions::Orientation orient,
1604  const int nq0,
1605  Array<OneD, int> &idmap)
1606  {
1607  boost::ignore_unused(nvert);
1608 
1609  if (idmap.num_elements() != nq0)
1610  {
1611  idmap = Array<OneD, int>(nq0);
1612  }
1613  switch (orient)
1614  {
1615  case StdRegions::eForwards:
1616  // Fwd
1617  for (int i = 0; i < nq0; ++i)
1618  {
1619  idmap[i] = i;
1620  }
1621  break;
1622  case StdRegions::eBackwards:
1623  {
1624  // Bwd
1625  for (int i = 0; i < nq0; ++i)
1626  {
1627  idmap[i] = nq0-1-i;
1628  }
1629  }
1630  break;
1631  default:
1632  ASSERTL0(false, "Unknown orientation");
1633  break;
1634  }
1635  }
1636 
1637  Array<OneD, NekDouble> Expansion2D::v_GetMF(
1638  const int dir,
1639  const int shapedim,
1640  const StdRegions::VarCoeffMap &varcoeffs)
1641  {
1642  return Expansion::v_GetMF(dir,shapedim,varcoeffs);
1643  }
1644 
1645  Array<OneD, NekDouble> Expansion2D::v_GetMFDiv(
1646  const int dir,
1647  const StdRegions::VarCoeffMap &varcoeffs)
1648  {
1649  return Expansion::v_GetMFDiv(dir,varcoeffs);
1650  }
1651 
1652  Array<OneD, NekDouble> Expansion2D::v_GetMFMag(
1653  const int dir,
1654  const StdRegions::VarCoeffMap &varcoeffs)
1655  {
1656  return Expansion::v_GetMFMag(dir,varcoeffs);
1657  }
1658 
1659  // Compute edgenormal \cdot vector
1660  Array<OneD, NekDouble> Expansion2D::v_GetnEdgecdotMF(
1661  const int dir,
1662  const int edge,
1663  ExpansionSharedPtr &EdgeExp_e,
1664  const Array<OneD, const Array<OneD, NekDouble> > &normals,
1665  const StdRegions::VarCoeffMap &varcoeffs)
1666  {
1667  int nquad_e = EdgeExp_e->GetNumPoints(0);
1668  int coordim = GetCoordim();
1669  int nquad0 = m_base[0]->GetNumPoints();
1670  int nquad1 = m_base[1]->GetNumPoints();
1671  int nqtot = nquad0*nquad1;
1672 
1673  StdRegions::VarCoeffType MMFCoeffs[15] =
1674  {
1675  StdRegions::eVarCoeffMF1x,
1676  StdRegions::eVarCoeffMF1y,
1677  StdRegions::eVarCoeffMF1z,
1678  StdRegions::eVarCoeffMF1Div,
1679  StdRegions::eVarCoeffMF1Mag,
1680  StdRegions::eVarCoeffMF2x,
1681  StdRegions::eVarCoeffMF2y,
1682  StdRegions::eVarCoeffMF2z,
1683  StdRegions::eVarCoeffMF2Div,
1684  StdRegions::eVarCoeffMF2Mag,
1685  StdRegions::eVarCoeffMF3x,
1686  StdRegions::eVarCoeffMF3y,
1687  StdRegions::eVarCoeffMF3z,
1688  StdRegions::eVarCoeffMF3Div,
1689  StdRegions::eVarCoeffMF3Mag
1690  };
1691 
1692  StdRegions::VarCoeffMap::const_iterator MFdir;
1693 
1694  Array<OneD, NekDouble> ncdotMF(nqtot,0.0);
1695  Array<OneD, NekDouble> tmp(nqtot);
1696  Array<OneD, NekDouble> tmp_e(nquad_e);
1697  for (int k=0; k<coordim; k++)
1698  {
1699  MFdir = varcoeffs.find(MMFCoeffs[dir*5+k]);
1700  tmp = MFdir->second;
1701 
1702  GetPhysEdgeVarCoeffsFromElement(edge, EdgeExp_e, tmp, tmp_e);
1703 
1704  Vmath::Vvtvp(nquad_e,
1705  &tmp_e[0], 1,
1706  &normals[k][0], 1,
1707  &ncdotMF[0], 1,
1708  &ncdotMF[0], 1);
1709  }
1710  return ncdotMF;
1711  }
1712 
1713  NekDouble Expansion2D::v_VectorFlux(
1714  const Array<OneD, Array<OneD, NekDouble> > &vec)
1715  {
1716  const Array<OneD, const Array<OneD, NekDouble> >
1717  &normals = GetLeftAdjacentElementExp()->
1718  GetFaceNormal(GetLeftAdjacentElementFace());
1719 
1720  int nq = GetTotPoints();
1721  Array<OneD, NekDouble > Fn(nq);
1722  Vmath::Vmul (nq, &vec[0][0], 1, &normals[0][0], 1, &Fn[0], 1);
1723  Vmath::Vvtvp(nq, &vec[1][0], 1, &normals[1][0], 1, &Fn[0], 1, &Fn[0], 1);
1724  Vmath::Vvtvp(nq, &vec[2][0], 1, &normals[2][0], 1, &Fn[0], 1, &Fn[0], 1);
1725 
1726  return StdExpansion::Integral(Fn);
1727  }
1728  }
1729 }
Nektar::StdRegions::StdMatrixKey::GetVarCoeffs
const VarCoeffMap & GetVarCoeffs() const
Definition: StdMatrixKey.h:161
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::LocalRegions::Expansion2D::v_GetMFMag
virtual Array< OneD, NekDouble > v_GetMFMag(const int dir, const StdRegions::VarCoeffMap &varcoeffs)
Definition: Expansion2D.cpp:1652
Nektar::LocalRegions::Expansion2D::GetLeftAdjacentElementExp
Expansion3DSharedPtr GetLeftAdjacentElementExp() const
Definition: Expansion2D.h:245
Nektar::SpatialDomains::Geometry2DSharedPtr
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry.h:65
Nektar::StdRegions::StdExpansion::DetShapeType
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:469
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:16
Nektar::DNekScalMatSharedPtr
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Definition: NekMatrixFwd.hpp:69
Nektar::LocalRegions::Expansion2D::GetPhysEdgeVarCoeffsFromElement
void GetPhysEdgeVarCoeffsFromElement(const int edge, ExpansionSharedPtr &EdgeExp, const Array< OneD, const NekDouble > &varcoeff, Array< OneD, NekDouble > &outarray)
Definition: Expansion2D.cpp:678
Nektar::LocalRegions::Expansion2D::AddNormTraceInt
void AddNormTraceInt(const int dir, Array< OneD, ExpansionSharedPtr > &EdgeExp, Array< OneD, Array< OneD, NekDouble > > &edgeCoeffs, Array< OneD, NekDouble > &outarray)
Definition: Expansion2D.cpp:402
Nektar::StdRegions::StdMatrixKey::GetVarCoeffAsMap
const VarCoeffMap GetVarCoeffAsMap(const VarCoeffType &coeff) const
Definition: StdMatrixKey.h:154
Nektar::StdRegions::StdExpansion::GetEdgeInteriorMap
void GetEdgeInteriorMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray)
Definition: StdExpansion.h:828
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:177
Nektar::LibUtilities::InterpCoeff1D
void InterpCoeff1D(const BasisKey &fbasis0, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, Array< OneD, NekDouble > &to)
Definition: InterpCoeff.cpp:46
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:634
Nektar::LocalRegions::Expansion2D::v_GetnEdgecdotMF
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)
Definition: Expansion2D.cpp:1660
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
Geometry.h
Nektar::LocalRegions::Expansion2D::m_negatedNormals
std::map< int, bool > m_negatedNormals
Definition: Expansion2D.h:135
Expansion2D.h
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:488
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:445
Nektar::LibUtilities::eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
Nektar::LocalRegions::Expansion2D::m_edgeNormals
std::map< int, StdRegions::NormalVector > m_edgeNormals
Definition: Expansion2D.h:134
Nektar::LibUtilities::eGauss_Lagrange
Lagrange Polynomials using the Gauss points .
Definition: BasisType.h:55
Nektar::StdRegions::StdExpansion::GetEorient
StdRegions::Orientation GetEorient(int edge)
Definition: StdExpansion.h:776
std
STL namespace.
Nektar::LocalRegions::Expansion2D::SetTraceToGeomOrientation
void SetTraceToGeomOrientation(Array< OneD, ExpansionSharedPtr > &EdgeExp, Array< OneD, NekDouble > &inout)
Definition: Expansion2D.cpp:309
Nektar::LocalRegions::Expansion2D::v_AddRobinEdgeContribution
virtual void v_AddRobinEdgeContribution(const int edgeid, const Array< OneD, const NekDouble > &primCoeffs, Array< OneD, NekDouble > &coeffs)
Definition: Expansion2D.cpp:1461
Nektar::StdRegions::ConstFactorMap
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:294
Nektar::LocalRegions::Expansion2D::v_BuildVertexMatrix
virtual DNekMatSharedPtr v_BuildVertexMatrix(const DNekScalMatSharedPtr &r_bnd)
Definition: Expansion2D.cpp:1496
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:69
Nektar::StdRegions::StdMatrixKey::GetMatrixType
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:81
Nektar::StdRegions::StdExpansion::GetNverts
int GetNverts() const
This function returns the number of vertices of the expansion domain.
Definition: StdExpansion.h:258
Nektar::StdRegions::StdMatrixKey::GetVarCoeff
const Array< OneD, const NekDouble > & GetVarCoeff(const StdRegions::VarCoeffType &coeff) const
Definition: StdMatrixKey.h:145
Nektar::LocalRegions::Expansion2D::v_AddRobinMassMatrix
virtual void v_AddRobinMassMatrix(const int edgeid, const Array< OneD, const NekDouble > &primCoeffs, DNekMatSharedPtr &inoutmat)
Definition: Expansion2D.cpp:1336
Nektar::StdRegions::StdMatrixKey::HasVarCoeff
bool HasVarCoeff(const StdRegions::VarCoeffType &coeff) const
Definition: StdMatrixKey.h:166
Nektar::LocalRegions::Expansion2D::AddEdgeBoundaryInt
void AddEdgeBoundaryInt(const int edge, ExpansionSharedPtr &EdgeExp, Array< OneD, NekDouble > &edgePhys, Array< OneD, NekDouble > &outarray, const StdRegions::VarCoeffMap &varcoeffs=StdRegions::NullVarCoeffMap)
Definition: Expansion2D.cpp:437
Geometry2D.h
Nektar::StdRegions::StdExpansion::GetFaceNormal
const NormalVector & GetFaceNormal(const int face) const
Definition: StdExpansion.h:1315
Nektar::LocalRegions::Expansion2D::v_NegateEdgeNormal
virtual void v_NegateEdgeNormal(const int edge)
Definition: Expansion2D.cpp:1586
Nektar::LocalRegions::Expansion2D::m_requireNeg
std::vector< bool > m_requireNeg
Definition: Expansion2D.h:133
Nektar::LocalRegions::Expansion2D::v_GetMF
virtual Array< OneD, NekDouble > v_GetMF(const int dir, const int shapedim, const StdRegions::VarCoeffMap &varcoeffs)
Definition: Expansion2D.cpp:1637
Nektar::StdRegions::StdExpansion::GetEdgeNcoeffs
int GetEdgeNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th edge.
Definition: StdExpansion.h:286
Nektar::StdRegions::StdExpansion::GetVertexMap
int GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
Definition: StdExpansion.h:822
Nektar::LocalRegions::Expansion2D::GetGeom2D
SpatialDomains::Geometry2DSharedPtr GetGeom2D() const
Definition: Expansion2D.h:291
Nektar::StdRegions::VarCoeffMap
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
Definition: StdRegions.hpp:264
ManagerAccess.h
Nektar::LocalRegions::Expansion2D::v_GetMFDiv
virtual Array< OneD, NekDouble > v_GetMFDiv(const int dir, const StdRegions::VarCoeffMap &varcoeffs)
Definition: Expansion2D.cpp:1645
Nektar::LocalRegions::Expansion::v_GetMFMag
Array< OneD, NekDouble > v_GetMFMag(const int dir, const StdRegions::VarCoeffMap &varcoeffs)
Definition: Expansion.cpp:444
Nektar::LocalRegions::Expansion2D::GetEdgeExp
Expansion1DSharedPtr GetEdgeExp(int edge, bool SetUpNormal=true)
Definition: Expansion2D.h:222
Nektar::LocalRegions::Expansion::v_GetMF
Array< OneD, NekDouble > v_GetMF(const int dir, const int shapedim, const StdRegions::VarCoeffMap &varcoeffs)
Definition: Expansion.cpp:347
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:161
Nektar::Transpose
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Definition: BlockMatrix.hpp:217
Expansion3D.h
Nektar::LocalRegions::Expansion2D::v_AddEdgeNormBoundaryInt
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)
Definition: Expansion2D.cpp:59
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:399
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::StdRegions::StdExpansion::GetNedges
int GetNedges() const
This function returns the number of edges of the expansion domain.
Definition: StdExpansion.h:271
Nektar::LocalRegions::Expansion2D::AddHDGHelmholtzEdgeTerms
void AddHDGHelmholtzEdgeTerms(const NekDouble tau, const int edge, Array< OneD, ExpansionSharedPtr > &EdgeExp, Array< OneD, NekDouble > &edgePhys, const StdRegions::VarCoeffMap &dirForcing, Array< OneD, NekDouble > &outarray)
Definition: Expansion2D.cpp:511
Nektar::LocalRegions::Expansion2D::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: Expansion2D.cpp:713
Nektar::LocalRegions::Expansion2D::v_EdgeNormalNegated
virtual bool v_EdgeNormalNegated(const int edge)
Definition: Expansion2D.cpp:1596
Nektar::StdRegions::StdMatrixKey::GetConstFactor
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:121
MatrixKey.h
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::LocalRegions::ExpansionSharedPtr
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:65
Nektar::LocalRegions::Expansion::v_GetMFDiv
Array< OneD, NekDouble > v_GetMFDiv(const int dir, const StdRegions::VarCoeffMap &varcoeffs)
Definition: Expansion.cpp:411
Nektar::LocalRegions::Expansion::GetLocMatrix
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition: Expansion.cpp:85
Nektar::LocalRegions::Expansion2D::AddHDGHelmholtzTraceTerms
void AddHDGHelmholtzTraceTerms(const NekDouble tau, const Array< OneD, const NekDouble > &inarray, Array< OneD, ExpansionSharedPtr > &EdgeExp, const StdRegions::VarCoeffMap &dirForcing, Array< OneD, NekDouble > &outarray)
Definition: Expansion2D.cpp:479
Nektar::StdRegions::StdExpansion::GetEdgeNormal
const NormalVector & GetEdgeNormal(const int edge) const
Definition: StdExpansion.h:1265
Nektar::StdRegions::StdMatrixKey::GetConstFactors
const ConstFactorMap & GetConstFactors() const
Definition: StdMatrixKey.h:135
Nektar::LocalRegions::Expansion::AddEdgeNormBoundaryInt
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:115
Nektar::StdRegions::StdExpansion::GetEdgeToElementMap
void GetEdgeToElementMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, int P=-1)
Definition: StdExpansion.h:849
Nektar::LocalRegions::Expansion2D::v_DGDeriv
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)
Definition: Expansion2D.cpp:1298
Nektar::LocalRegions::Expansion2D::m_edgeExp
std::vector< Expansion1DWeakPtr > m_edgeExp
Definition: Expansion2D.h:132
Nektar::StdRegions::StdExpansion::v_GetEorient
virtual StdRegions::Orientation v_GetEorient(int edge)
Definition: StdExpansion.cpp:1060
Nektar::LocalRegions::Expansion2D::v_GetSurfaceNormal
const StdRegions::NormalVector & v_GetSurfaceNormal(const int id) const
Definition: Expansion2D.cpp:1580
Nektar::LocalRegions::Expansion2D::v_GetEdgeNormal
const StdRegions::NormalVector & v_GetEdgeNormal(const int edge) const
Definition: Expansion2D.cpp:1571
Nektar::LocalRegions::Expansion2D::v_SetUpPhysNormals
virtual void v_SetUpPhysNormals(const int edge)
Definition: Expansion2D.cpp:1566
Nektar::StdRegions::StdExpansion::m_stdMatrixManager
LibUtilities::NekManager< StdMatrixKey, DNekMat, StdMatrixKey::opLess > m_stdMatrixManager
Definition: StdExpansion.h:1432
Nektar::StdRegions::StdExpansion::GetTotPoints
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:140
Expansion1D.h
InterpCoeff.h
Nektar::StdRegions::IndexMapValuesSharedPtr
std::shared_ptr< IndexMapValues > IndexMapValuesSharedPtr
Definition: IndexMapKey.h:124
Nektar::LibUtilities::eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:54
Nektar::LocalRegions::Expansion2D::v_VectorFlux
virtual NekDouble v_VectorFlux(const Array< OneD, Array< OneD, NekDouble > > &vec)
Definition: Expansion2D.cpp:1713
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:376
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:250
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition: StdExpansion.h:1429
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
Nektar::LocalRegions::Expansion1DSharedPtr
std::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:51
Nektar::StdRegions::StdExpansion::GetBoundaryMap
void GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:812
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:49
Nektar::LocalRegions::Expansion2D::v_GetEdgeInverseBoundaryMap
Array< OneD, unsigned int > v_GetEdgeInverseBoundaryMap(int eid)
Definition: Expansion2D.cpp:1527
Nektar::LocalRegions::Expansion2D::ReOrientEdgePhysMap
void ReOrientEdgePhysMap(const int nvert, const StdRegions::Orientation orient, const int nq0, Array< OneD, int > &idmap)
Definition: Expansion2D.cpp:1601
Nektar::StdRegions::StdExpansion::ComputeEdgeNormal
void ComputeEdgeNormal(const int edge)
Definition: StdExpansion.h:1270
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:186
Nektar::StdRegions::NullConstFactorMap
static ConstFactorMap NullConstFactorMap
Definition: StdRegions.hpp:295
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:1957
Generated on Wed Jan 20 2021 03:33:33 for Nektar++ by   doxygen 1.8.13