Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: File for Expansion2D routines
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 #include <LocalRegions/Expansion2D.h>
36 #include <LocalRegions/Expansion1D.h>
37 #include <SpatialDomains/Geometry.h>
38 #include <SpatialDomains/Geometry2D.h>
39 #include <LibUtilities/Foundations/InterpCoeff.h>
40 #include <LocalRegions/MatrixKey.h>
41 #include <LibUtilities/Foundations/ManagerAccess.h>
42 
43 namespace Nektar
44 {
45  namespace LocalRegions
46  {
47  Expansion2D::Expansion2D(SpatialDomains::Geometry2DSharedPtr pGeom) :
48  StdExpansion(), Expansion(pGeom), StdExpansion2D()
49  {
50  m_elementFaceLeft = -1;
51  m_elementFaceRight = -1;
52  }
53 
54  void Expansion2D::v_AddEdgeNormBoundaryInt(
55 
56  const int edge,
57 
58  StdRegions::StdExpansionSharedPtr &EdgeExp,
59  const Array<OneD, const NekDouble> &Fx,
60  const Array<OneD, const NekDouble> &Fy,
61  Array<OneD, NekDouble> &outarray)
62  {
63  ASSERTL1(GetCoordim() == 2,
64  "Routine only set up for two-dimensions");
65 
66  const Array<OneD, const Array<OneD, NekDouble> > normals
67  = GetEdgeNormal(edge);
68 
69  if (m_requireNeg.size() == 0)
70  {
71  m_requireNeg.resize(GetNedges());
72 
73  for (int i = 0; i < GetNedges(); ++i)
74  {
75  m_requireNeg[i] = false;
76  if (m_negatedNormals[i])
77  {
78  m_requireNeg[i] = true;
79  continue;
80  }
81 
82  Expansion1DSharedPtr edgeExp =
83  m_edgeExp[i].lock()->as<Expansion1D>();
84 
85  if (edgeExp->GetRightAdjacentElementExp())
86  {
87  if (edgeExp->GetRightAdjacentElementExp()->GetGeom2D()
88  ->GetGlobalID() == GetGeom2D()->GetGlobalID())
89  {
90  m_requireNeg[i] = true;
91  }
92  }
93  }
94  }
95 
96  // We allow the case of mixed polynomial order by supporting only
97  // those modes on the edge common to both adjoining elements. This
98  // is enforced here by taking the minimum size and padding with
99  // zeros.
100  int nquad_e = min(EdgeExp->GetNumPoints(0),
101  int(normals[0].num_elements()));
102 
103  int nEdgePts = EdgeExp->GetTotPoints();
104  Array<OneD, NekDouble> edgePhys(nEdgePts);
105  Vmath::Vmul (nquad_e, normals[0], 1, Fx, 1, edgePhys, 1);
106  Vmath::Vvtvp(nquad_e, normals[1], 1, Fy, 1, edgePhys, 1,
107  edgePhys, 1);
108 
109  Expansion1DSharedPtr locExp = EdgeExp->as<Expansion1D>();
110 
111  if (m_negatedNormals[edge])
112  {
113  Vmath::Neg(nquad_e, edgePhys, 1);
114  }
115  else if (locExp->GetRightAdjacentElementEdge() != -1)
116  {
117  if (locExp->GetRightAdjacentElementExp()->GetGeom2D()->GetGlobalID()
118  == GetGeom2D()->GetGlobalID())
119  {
120  Vmath::Neg(nquad_e, edgePhys, 1);
121  }
122  }
123 
124  AddEdgeNormBoundaryInt(edge, EdgeExp, edgePhys, outarray);
125  }
126 
127  void Expansion2D::v_AddEdgeNormBoundaryInt(
128  const int edge,
129  StdRegions::StdExpansionSharedPtr &EdgeExp,
130  const Array<OneD, const NekDouble> &Fn,
131  Array<OneD, NekDouble> &outarray)
132  {
133  int i;
134 
135  if (m_requireNeg.size() == 0)
136  {
137  m_requireNeg.resize(GetNedges());
138 
139  for (i = 0; i < GetNedges(); ++i)
140  {
141  m_requireNeg[i] = false;
142  if (m_negatedNormals[i])
143  {
144  m_requireNeg[i] = true;
145  continue;
146  }
147 
148  Expansion1DSharedPtr edgeExp =
149  m_edgeExp[i].lock()->as<Expansion1D>();
150 
151  if (edgeExp->GetRightAdjacentElementExp())
152  {
153  if (edgeExp->GetRightAdjacentElementExp()->GetGeom2D()
154  ->GetGlobalID() == GetGeom2D()->GetGlobalID())
155  {
156  m_requireNeg[i] = true;
157  }
158  }
159  }
160  }
161 
162  StdRegions::IndexMapKey ikey(
163  StdRegions::eEdgeToElement, DetShapeType(),
164  GetBasisNumModes(0), GetBasisNumModes(1), 0,
165  edge, GetEorient(edge));
166  StdRegions::IndexMapValuesSharedPtr map =
167  StdExpansion::GetIndexMap(ikey);
168 
169  // Order of the element
170  int order_e = map->num_elements();
171  // Order of the trace
172  int n_coeffs = EdgeExp->GetNcoeffs();
173 
174  Array<OneD, NekDouble> edgeCoeffs(n_coeffs);
175  if(n_coeffs!=order_e) // Going to orthogonal space
176  {
177  EdgeExp->FwdTrans(Fn, edgeCoeffs);
178  Expansion1DSharedPtr locExp = EdgeExp->as<Expansion1D>();
179 
180  if (m_requireNeg[edge])
181  {
182  Vmath::Neg(n_coeffs, edgeCoeffs, 1);
183  }
184 
185  Array<OneD, NekDouble> coeff(n_coeffs,0.0);
186  LibUtilities::BasisType btype = ((LibUtilities::BasisType) 1); //1-->Ortho_A
187  LibUtilities::BasisKey bkey_ortho(btype,EdgeExp->GetBasis(0)->GetNumModes(),EdgeExp->GetBasis(0)->GetPointsKey());
188  LibUtilities::BasisKey bkey(EdgeExp->GetBasis(0)->GetBasisType(),EdgeExp->GetBasis(0)->GetNumModes(),EdgeExp->GetBasis(0)->GetPointsKey());
189  LibUtilities::InterpCoeff1D(bkey,edgeCoeffs,bkey_ortho,coeff);
190 
191  // Cutting high frequencies
192  for(i = order_e; i < n_coeffs; i++)
193  {
194  coeff[i] = 0.0;
195  }
196 
197  LibUtilities::InterpCoeff1D(bkey_ortho,coeff,bkey,edgeCoeffs);
198 
199  StdRegions::StdMatrixKey masskey(StdRegions::eMass,LibUtilities::eSegment,*EdgeExp);
200  EdgeExp->MassMatrixOp(edgeCoeffs, edgeCoeffs, masskey);
201  }
202  else
203  {
204  EdgeExp->IProductWRTBase(Fn, edgeCoeffs);
205 
206  Expansion1DSharedPtr locExp = EdgeExp->as<Expansion1D>();
207 
208  if (m_requireNeg[edge])
209  {
210  Vmath::Neg(n_coeffs, edgeCoeffs, 1);
211  }
212  }
213 
214  // Implementation for all the basis except Gauss points
215  if(EdgeExp->GetBasis(0)->GetBasisType() !=
216  LibUtilities::eGauss_Lagrange)
217  {
218  // add data to outarray if forward edge normal is outwards
219  for(i = 0; i < order_e; ++i)
220  {
221  outarray[(*map)[i].index] +=
222  (*map)[i].sign * edgeCoeffs[i];
223  }
224  }
225  else
226  {
227  int nCoeffs0, nCoeffs1;
228  int j;
229 
230  StdRegions::ConstFactorMap factors;
231  factors[StdRegions::eFactorGaussEdge] = edge;
232  StdRegions::StdMatrixKey key(StdRegions::eGaussDG,
233  DetShapeType(),*this,factors);
234 
235  DNekMatSharedPtr mat_gauss = m_stdMatrixManager[key];
236 
237  switch(edge)
238  {
239  case 0:
240  {
241  nCoeffs1 = m_base[1]->GetNumModes();
242 
243  for(i = 0; i < order_e; ++i)
244  {
245  for(j = 0; j < nCoeffs1; j++)
246  {
247  outarray[(*map)[i].index + j*order_e] +=
248  mat_gauss->GetPtr()[j]*
249  (*map)[i].sign*edgeCoeffs[i];
250  }
251  }
252  break;
253  }
254  case 1:
255  {
256  nCoeffs0 = m_base[0]->GetNumModes();
257 
258  for(i = 0; i < order_e; ++i)
259  {
260  for(j = 0; j < nCoeffs0; j++)
261  {
262  outarray[(*map)[i].index - j] +=
263  mat_gauss->GetPtr()[order_e - 1 -j]*
264  (*map)[i].sign*edgeCoeffs[i];
265  }
266  }
267  break;
268  }
269  case 2:
270  {
271  nCoeffs1 = m_base[1]->GetNumModes();
272 
273  for(i = 0; i < order_e; ++i)
274  {
275  for(j = 0; j < nCoeffs1; j++)
276  {
277  outarray[(*map)[i].index - j*order_e] +=
278  mat_gauss->GetPtr()[order_e - 1 - j]*
279  (*map)[i].sign*edgeCoeffs[i];
280  }
281  }
282  break;
283  }
284  case 3:
285  {
286  nCoeffs0 = m_base[0]->GetNumModes();
287 
288  for(i = 0; i < order_e; ++i)
289  {
290  for(j = 0; j < nCoeffs0; j++)
291  {
292  outarray[(*map)[i].index + j] +=
293  mat_gauss->GetPtr()[j]*
294  (*map)[i].sign*edgeCoeffs[i];
295  }
296  }
297  break;
298  }
299  default:
300  ASSERTL0(false,"edge value (< 3) is out of range");
301  break;
302  }
303  }
304  }
305 
306  void Expansion2D::SetTraceToGeomOrientation(
307  Array<OneD,StdRegions::StdExpansionSharedPtr> &EdgeExp,
308  Array<OneD, NekDouble> &inout)
309  {
310  int i, cnt = 0;
311  int nedges = GetNedges();
312  Array<OneD, NekDouble> e_tmp;
313 
314  for(i = 0; i < nedges; ++i)
315  {
316  EdgeExp[i]->SetCoeffsToOrientation(GetEorient(i),
317  e_tmp = inout + cnt,
318  e_tmp = inout + cnt);
319  cnt += GetEdgeNcoeffs(i);
320  }
321  }
322 
323  /**
324  * Computes the C matrix entries due to the presence of the identity
325  * matrix in Eqn. 32.
326  */
327  void Expansion2D::AddNormTraceInt(
328  const int dir,
329  Array<OneD, const NekDouble> &inarray,
330  Array<OneD, StdRegions::StdExpansionSharedPtr> &EdgeExp,
331  Array<OneD, NekDouble> &outarray,
332  const StdRegions::VarCoeffMap &varcoeffs)
333  {
334  int i,e,cnt;
335  int order_e,nquad_e;
336  int nedges = GetNedges();
337 
338  cnt = 0;
339  for(e = 0; e < nedges; ++e)
340  {
341  order_e = EdgeExp[e]->GetNcoeffs();
342  nquad_e = EdgeExp[e]->GetNumPoints(0);
343 
344  const Array<OneD, const Array<OneD, NekDouble> > &normals
345  = GetEdgeNormal(e);
346  Array<OneD, NekDouble> edgeCoeffs(order_e);
347  Array<OneD, NekDouble> edgePhys (nquad_e);
348 
349  for(i = 0; i < order_e; ++i)
350  {
351  edgeCoeffs[i] = inarray[i+cnt];
352  }
353  cnt += order_e;
354 
355  EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
356 
357  // Multiply by variable coefficient
358  /// @TODO: Document this
359  // StdRegions::VarCoeffType VarCoeff[3] = {StdRegions::eVarCoeffD00,
360  // StdRegions::eVarCoeffD11,
361  // StdRegions::eVarCoeffD22};
362  // StdRegions::VarCoeffMap::const_iterator x;
363  // Array<OneD, NekDouble> varcoeff_work(nquad_e);
364 
365  // if ((x = varcoeffs.find(VarCoeff[dir])) != varcoeffs.end())
366  // {
367  // GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
368  // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
369  // }
370 
371  Vmath::Vmul(nquad_e, normals[dir], 1, edgePhys, 1, edgePhys, 1);
372 
373  if (m_negatedNormals[e])
374  {
375  Vmath::Neg(nquad_e, edgePhys, 1);
376  }
377 
378  AddEdgeBoundaryInt(e, EdgeExp[e], edgePhys, outarray, varcoeffs);
379  }
380  }
381 
382  void Expansion2D::AddNormTraceInt(
383  const int dir,
384  Array<OneD, StdRegions::StdExpansionSharedPtr> &EdgeExp,
385  Array<OneD, Array<OneD, NekDouble> > &edgeCoeffs,
386  Array<OneD, NekDouble> &outarray)
387  {
388  int e;
389  int nquad_e;
390  int nedges = GetNedges();
391 
392  for(e = 0; e < nedges; ++e)
393  {
394  nquad_e = EdgeExp[e]->GetNumPoints(0);
395 
396  Array<OneD, NekDouble> edgePhys(nquad_e);
397  const Array<OneD, const Array<OneD, NekDouble> > &normals
398  = GetEdgeNormal(e);
399 
400  EdgeExp[e]->BwdTrans(edgeCoeffs[e], edgePhys);
401 
402  Vmath::Vmul(nquad_e, normals[dir], 1, edgePhys, 1, edgePhys, 1);
403 
404  if (m_negatedNormals[e])
405  {
406  Vmath::Neg(nquad_e, edgePhys, 1);
407  }
408 
409  AddEdgeBoundaryInt(e, EdgeExp[e], edgePhys, outarray);
410  }
411  }
412 
413 
414  /**
415  * For a given edge add the \tilde{F}_1j contributions
416  */
417  void Expansion2D::AddEdgeBoundaryInt(
418  const int edge,
419  StdRegions::StdExpansionSharedPtr &EdgeExp,
420  Array<OneD, NekDouble> &edgePhys,
421  Array<OneD, NekDouble> &outarray,
422  const StdRegions::VarCoeffMap &varcoeffs)
423  {
424  int i;
425  int order_e = EdgeExp->GetNcoeffs();
426  int nquad_e = EdgeExp->GetNumPoints(0);
427  Array<OneD,unsigned int> map;
428  Array<OneD,int> sign;
429  Array<OneD, NekDouble> coeff(order_e);
430 
431  GetEdgeToElementMap(edge, v_GetEorient(edge), map, sign);
432 
433  StdRegions::VarCoeffType VarCoeff[3] = {StdRegions::eVarCoeffD00,
434  StdRegions::eVarCoeffD11,
435  StdRegions::eVarCoeffD22};
436  StdRegions::VarCoeffMap::const_iterator x;
437 
438  /// @TODO Variable coeffs
439  if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
440  {
441  Array<OneD, NekDouble> work(nquad_e);
442  GetPhysEdgeVarCoeffsFromElement(
443  edge, EdgeExp, x->second, work);
444  Vmath::Vmul(nquad_e, work, 1, edgePhys, 1, edgePhys, 1);
445  }
446 
447  EdgeExp->IProductWRTBase(edgePhys, coeff);
448 
449  // add data to out array
450  for(i = 0; i < order_e; ++i)
451  {
452  outarray[map[i]] += sign[i]*coeff[i];
453  }
454  }
455 
456  // This method assumes that data in EdgeExp is orientated according to
457  // elemental counter clockwise format AddHDGHelmholtzTraceTerms with
458  // directions
459  void Expansion2D::AddHDGHelmholtzTraceTerms(
460  const NekDouble tau,
461  const Array<OneD, const NekDouble> &inarray,
462  Array<OneD,StdRegions::StdExpansionSharedPtr> &EdgeExp,
463  const StdRegions::VarCoeffMap &dirForcing,
464  Array<OneD,NekDouble> &outarray)
465  {
466  ASSERTL0(&inarray[0] != &outarray[0],
467  "Input and output arrays use the same memory");
468 
469  int e, cnt, order_e, nedges = GetNedges();
470  Array<OneD, const NekDouble> tmp;
471 
472  cnt = 0;
473 
474  for(e = 0; e < nedges; ++e)
475  {
476  order_e = EdgeExp[e]->GetNcoeffs();
477  Array<OneD, NekDouble> edgeCoeffs(order_e);
478  Array<OneD, NekDouble> edgePhys (EdgeExp[e]->GetTotPoints());
479 
480  Vmath::Vcopy(order_e, tmp = inarray + cnt, 1, edgeCoeffs, 1);
481  EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
482  AddHDGHelmholtzEdgeTerms(
483  tau, e, EdgeExp, edgePhys, dirForcing, outarray);
484 
485  cnt += order_e;
486  }
487  }
488 
489  // evaluate additional terms in HDG edges. Not that this assumes that
490  // edges are unpacked into local cartesian order.
491  void Expansion2D::AddHDGHelmholtzEdgeTerms(
492  const NekDouble tau,
493  const int edge,
494  Array<OneD, StdRegions::StdExpansionSharedPtr> &EdgeExp,
495  Array<OneD, NekDouble> &edgePhys,
496  const StdRegions::VarCoeffMap &varcoeffs,
497  Array<OneD, NekDouble> &outarray)
498  {
499  int i, j, n;
500  int nquad_e = EdgeExp[edge]->GetNumPoints(0);
501  int order_e = EdgeExp[edge]->GetNcoeffs();
502  int coordim = GetCoordim();
503  int ncoeffs = GetNcoeffs();
504 
505  Array<OneD, NekDouble> inval (nquad_e);
506  Array<OneD, NekDouble> outcoeff(order_e);
507  Array<OneD, NekDouble> tmpcoeff(ncoeffs);
508 
509  const Array<OneD, const Array<OneD, NekDouble> > &normals
510  = GetEdgeNormal(edge);
511 
512  Array<OneD,unsigned int> emap;
513  Array<OneD,int> sign;
514 
515  DNekScalMat &invMass = *GetLocMatrix(StdRegions::eInvMass);
516 
517  StdRegions::Orientation edgedir = GetEorient(edge);
518 
519  DNekVec Coeffs (ncoeffs,outarray,eWrapper);
520  DNekVec Tmpcoeff(ncoeffs,tmpcoeff,eWrapper);
521 
522  GetEdgeToElementMap(edge,edgedir,emap,sign);
523 
524  StdRegions::MatrixType DerivType[3] = {StdRegions::eWeakDeriv0,
525  StdRegions::eWeakDeriv1,
526  StdRegions::eWeakDeriv2};
527 
528  StdRegions::VarCoeffType VarCoeff[3] = {StdRegions::eVarCoeffD00,
529  StdRegions::eVarCoeffD11,
530  StdRegions::eVarCoeffD22};
531 
532  StdRegions::VarCoeffMap::const_iterator x;
533  /// @TODO: What direction to use here??
534  if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
535  {
536  Array<OneD, NekDouble> work(nquad_e);
537  GetPhysEdgeVarCoeffsFromElement(
538  edge, EdgeExp[edge], x->second, work);
539  Vmath::Vmul(nquad_e, work, 1, edgePhys, 1, edgePhys, 1);
540  }
541 
542  //================================================================
543  // Add F = \tau <phi_i,in_phys>
544  // Fill edge and take inner product
545  EdgeExp[edge]->IProductWRTBase(edgePhys, outcoeff);
546  // add data to out array
547  for(i = 0; i < order_e; ++i)
548  {
549  outarray[emap[i]] += sign[i] * tau * outcoeff[i];
550  }
551  //================================================================
552 
553  //===============================================================
554  // Add -\sum_i D_i^T M^{-1} G_i + E_i M^{-1} G_i =
555  // \sum_i D_i M^{-1} G_i term
556 
557  // Two independent direction
558  for(n = 0; n < coordim; ++n)
559  {
560  Vmath::Vmul(nquad_e, normals[n], 1, edgePhys, 1, inval, 1);
561 
562  if (m_negatedNormals[edge])
563  {
564  Vmath::Neg(nquad_e, inval, 1);
565  }
566 
567  // Multiply by variable coefficient
568  /// @TODO: Document this (probably not needed)
569 // StdRegions::VarCoeffMap::const_iterator x;
570 // if ((x = varcoeffs.find(VarCoeff[n])) != varcoeffs.end())
571 // {
572 // GetPhysEdgeVarCoeffsFromElement(edge,EdgeExp[edge],x->second,varcoeff_work);
573 // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[edge]->GetPhys(),1,EdgeExp[edge]->UpdatePhys(),1);
574 // }
575 
576  EdgeExp[edge]->IProductWRTBase(inval, outcoeff);
577 
578  // M^{-1} G
579  for(i = 0; i < ncoeffs; ++i)
580  {
581  tmpcoeff[i] = 0;
582  for(j = 0; j < order_e; ++j)
583  {
584  tmpcoeff[i] += invMass(i,emap[j])*sign[j]*outcoeff[j];
585  }
586  }
587 
588  if(varcoeffs.find(VarCoeff[n]) != varcoeffs.end())
589  {
590  MatrixKey mkey(DerivType[n], DetShapeType(), *this, StdRegions::NullConstFactorMap, varcoeffs);
591  DNekScalMat &Dmat = *GetLocMatrix(mkey);
592  Coeffs = Coeffs + Dmat*Tmpcoeff;
593  }
594 
595  else
596  {
597  DNekScalMat &Dmat = *GetLocMatrix(DerivType[n]);
598  Coeffs = Coeffs + Dmat*Tmpcoeff;
599  }
600  }
601  }
602 
603  /**
604  * Extracts the variable coefficients along an edge
605  */
606  void Expansion2D::GetPhysEdgeVarCoeffsFromElement(
607  const int edge,
608  StdRegions::StdExpansionSharedPtr &EdgeExp,
609  const Array<OneD, const NekDouble> &varcoeff,
610  Array<OneD,NekDouble> &outarray)
611  {
612  Array<OneD, NekDouble> tmp(GetNcoeffs());
613  Array<OneD, NekDouble> edgetmp(EdgeExp->GetNcoeffs());
614 
615  // FwdTrans varcoeffs
616  FwdTrans(varcoeff, tmp);
617 
618  // Map to edge
619  Array<OneD,unsigned int> emap;
620  Array<OneD, int> sign;
621  StdRegions::Orientation edgedir = GetEorient(edge);
622  GetEdgeToElementMap(edge,edgedir,emap,sign);
623 
624  for (unsigned int i = 0; i < EdgeExp->GetNcoeffs(); ++i)
625  {
626  edgetmp[i] = tmp[emap[i]];
627  }
628 
629  // BwdTrans
630  EdgeExp->BwdTrans(edgetmp, outarray);
631  }
632 
633 
634  /**
635  * Computes matrices needed for the HDG formulation. References to
636  * equations relate to the following paper:
637  * R. M. Kirby, S. J. Sherwin, B. Cockburn, To CG or to HDG: A
638  * Comparative Study, J. Sci. Comp P1-30
639  * DOI 10.1007/s10915-011-9501-7
640  */
641  DNekMatSharedPtr Expansion2D::v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
642  {
643  DNekMatSharedPtr returnval;
644 
645  switch(mkey.GetMatrixType())
646  {
647  // (Z^e)^{-1} (Eqn. 33, P22)
648  case StdRegions::eHybridDGHelmholtz:
649  {
650  ASSERTL1(IsBoundaryInteriorExpansion(),
651  "HybridDGHelmholtz matrix not set up "
652  "for non boundary-interior expansions");
653 
654  int i,j,k;
655  NekDouble lambdaval = mkey.GetConstFactor(StdRegions::eFactorLambda);
656  NekDouble tau = mkey.GetConstFactor(StdRegions::eFactorTau);
657  int ncoeffs = GetNcoeffs();
658  int nedges = GetNedges();
659 
660  Array<OneD,unsigned int> emap;
661  Array<OneD,int> sign;
662  StdRegions::Orientation edgedir = StdRegions::eForwards;
663  ExpansionSharedPtr EdgeExp;
664  StdRegions::StdExpansionSharedPtr EdgeExp2;
665 
666  int order_e, coordim = GetCoordim();
667  DNekScalMat &invMass = *GetLocMatrix(StdRegions::eInvMass);
668  StdRegions::MatrixType DerivType[3] = {StdRegions::eWeakDeriv0,
669  StdRegions::eWeakDeriv1,
670  StdRegions::eWeakDeriv2};
671  DNekMat LocMat(ncoeffs,ncoeffs);
672 
673  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(ncoeffs,ncoeffs);
674  DNekMat &Mat = *returnval;
675  Vmath::Zero(ncoeffs*ncoeffs,Mat.GetPtr(),1);
676 
677  StdRegions::VarCoeffType Coeffs[3] = {StdRegions::eVarCoeffD00,
678  StdRegions::eVarCoeffD11,
679  StdRegions::eVarCoeffD22};
680 
681  for(i=0; i < coordim; ++i)
682  {
683  if(mkey.HasVarCoeff(Coeffs[i]))
684  {
685  MatrixKey DmatkeyL(DerivType[i], DetShapeType(), *this, StdRegions::NullConstFactorMap, mkey.GetVarCoeffAsMap(Coeffs[i]));
686  MatrixKey DmatkeyR(DerivType[i], DetShapeType(), *this);
687 
688  DNekScalMat &DmatL = *GetLocMatrix(DmatkeyL);
689  DNekScalMat &DmatR = *GetLocMatrix(DmatkeyR);
690  Mat = Mat + DmatL*invMass*Transpose(DmatR);
691  }
692  else
693  {
694  DNekScalMat &Dmat = *GetLocMatrix(DerivType[i]);
695  Mat = Mat + Dmat*invMass*Transpose(Dmat);
696  }
697 
698  }
699 
700  // Add Mass Matrix Contribution for Helmholtz problem
701  DNekScalMat &Mass = *GetLocMatrix(StdRegions::eMass);
702  Mat = Mat + lambdaval*Mass;
703 
704  // Add tau*E_l using elemental mass matrices on each edge
705  for(i = 0; i < nedges; ++i)
706  {
707  EdgeExp = GetEdgeExp(i);
708  EdgeExp2 = GetEdgeExp(i);
709  order_e = EdgeExp->GetNcoeffs();
710  int nq = EdgeExp->GetNumPoints(0);
711  GetEdgeToElementMap(i,edgedir,emap,sign);
712 
713  // @TODO: Document
714  StdRegions::VarCoeffMap edgeVarCoeffs;
715  if (mkey.HasVarCoeff(StdRegions::eVarCoeffD00))
716  {
717  Array<OneD, NekDouble> mu(nq);
718  GetPhysEdgeVarCoeffsFromElement(i, EdgeExp2, mkey.GetVarCoeff(StdRegions::eVarCoeffD00), mu);
719  edgeVarCoeffs[StdRegions::eVarCoeffMass] = mu;
720  }
721  DNekScalMat &eMass = *EdgeExp->GetLocMatrix(StdRegions::eMass, StdRegions::NullConstFactorMap, edgeVarCoeffs);
722  //DNekScalMat &eMass = *EdgeExp->GetLocMatrix(StdRegions::eMass);
723 
724  for(j = 0; j < order_e; ++j)
725  {
726  for(k = 0; k < order_e; ++k)
727  {
728  Mat(emap[j],emap[k]) = Mat(emap[j],emap[k]) + tau*sign[j]*sign[k]*eMass(j,k);
729  }
730  }
731  }
732  }
733  break;
734  // U^e (P22)
735  case StdRegions::eHybridDGLamToU:
736  {
737  int i,j,k;
738  int nbndry = NumDGBndryCoeffs();
739  int ncoeffs = GetNcoeffs();
740  int nedges = GetNedges();
741  NekDouble tau = mkey.GetConstFactor(StdRegions::eFactorTau);
742 
743  Array<OneD,NekDouble> lambda(nbndry);
744  DNekVec Lambda(nbndry,lambda,eWrapper);
745  Array<OneD,NekDouble> ulam(ncoeffs);
746  DNekVec Ulam(ncoeffs,ulam,eWrapper);
747  Array<OneD,NekDouble> f(ncoeffs);
748  DNekVec F(ncoeffs,f,eWrapper);
749 
750  Array<OneD,StdRegions::StdExpansionSharedPtr> EdgeExp(nedges);
751  // declare matrix space
752  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(ncoeffs,nbndry);
753  DNekMat &Umat = *returnval;
754 
755  // Z^e matrix
756  MatrixKey newkey(StdRegions::eInvHybridDGHelmholtz, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
757  DNekScalMat &invHmat = *GetLocMatrix(newkey);
758 
759  Array<OneD,unsigned int> emap;
760  Array<OneD,int> sign;
761 
762  for(i = 0; i < nedges; ++i)
763  {
764  EdgeExp[i] = GetEdgeExp(i);
765  }
766 
767  // for each degree of freedom of the lambda space
768  // calculate Umat entry
769  // Generate Lambda to U_lambda matrix
770  for(j = 0; j < nbndry; ++j)
771  {
772  // standard basis vectors e_j
773  Vmath::Zero(nbndry,&lambda[0],1);
774  Vmath::Zero(ncoeffs,&f[0],1);
775  lambda[j] = 1.0;
776 
777  SetTraceToGeomOrientation(EdgeExp,lambda);
778 
779  // Compute F = [I D_1 M^{-1} D_2 M^{-1}] C e_j
780  AddHDGHelmholtzTraceTerms(
781  tau, lambda, EdgeExp, mkey.GetVarCoeffs(), f);
782 
783  // Compute U^e_j
784  Ulam = invHmat*F; // generate Ulam from lambda
785 
786  // fill column of matrix
787  for(k = 0; k < ncoeffs; ++k)
788  {
789  Umat(k,j) = Ulam[k];
790  }
791  }
792  }
793  break;
794  // Q_0, Q_1, Q_2 matrices (P23)
795  // Each are a product of a row of Eqn 32 with the C matrix.
796  // Rather than explicitly computing all of Eqn 32, we note each
797  // row is almost a multiple of U^e, so use that as our starting
798  // point.
799  case StdRegions::eHybridDGLamToQ0:
800  case StdRegions::eHybridDGLamToQ1:
801  case StdRegions::eHybridDGLamToQ2:
802  {
803  int i,j,k,dir;
804  int nbndry = NumDGBndryCoeffs();
805  int ncoeffs = GetNcoeffs();
806  int nedges = GetNedges();
807 
808  Array<OneD,NekDouble> lambda(nbndry);
809  DNekVec Lambda(nbndry,lambda,eWrapper);
810  Array<OneD,StdRegions::StdExpansionSharedPtr> EdgeExp(nedges);
811 
812  Array<OneD,NekDouble> ulam(ncoeffs);
813  DNekVec Ulam(ncoeffs,ulam,eWrapper);
814  Array<OneD,NekDouble> f(ncoeffs);
815  DNekVec F(ncoeffs,f,eWrapper);
816 
817  // declare matrix space
818  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(ncoeffs,nbndry);
819  DNekMat &Qmat = *returnval;
820 
821  // Lambda to U matrix
822  MatrixKey lamToUkey(StdRegions::eHybridDGLamToU, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
823  DNekScalMat &lamToU = *GetLocMatrix(lamToUkey);
824 
825  // Inverse mass matrix
826  DNekScalMat &invMass = *GetLocMatrix(StdRegions::eInvMass);
827 
828  for(i = 0; i < nedges; ++i)
829  {
830  EdgeExp[i] = GetEdgeExp(i);
831  }
832 
833  //Weak Derivative matrix
834  DNekScalMatSharedPtr Dmat;
835  switch(mkey.GetMatrixType())
836  {
837  case StdRegions::eHybridDGLamToQ0:
838  dir = 0;
839  Dmat = GetLocMatrix(StdRegions::eWeakDeriv0);
840  break;
841  case StdRegions::eHybridDGLamToQ1:
842  dir = 1;
843  Dmat = GetLocMatrix(StdRegions::eWeakDeriv1);
844  break;
845  case StdRegions::eHybridDGLamToQ2:
846  dir = 2;
847  Dmat = GetLocMatrix(StdRegions::eWeakDeriv2);
848  break;
849  default:
850  ASSERTL0(false,"Direction not known");
851  break;
852  }
853 
854  // for each degree of freedom of the lambda space
855  // calculate Qmat entry
856  // Generate Lambda to Q_lambda matrix
857  for(j = 0; j < nbndry; ++j)
858  {
859  Vmath::Zero(nbndry,&lambda[0],1);
860  lambda[j] = 1.0;
861 
862  // for lambda[j] = 1 this is the solution to ulam
863  for(k = 0; k < ncoeffs; ++k)
864  {
865  Ulam[k] = lamToU(k,j);
866  }
867 
868  // -D^T ulam
869  Vmath::Neg(ncoeffs,&ulam[0],1);
870  F = Transpose(*Dmat)*Ulam;
871 
872  SetTraceToGeomOrientation(EdgeExp,lambda);
873 
874  // Add the C terms resulting from the I's on the
875  // diagonals of Eqn 32
876  AddNormTraceInt(dir,lambda,EdgeExp,f,mkey.GetVarCoeffs());
877 
878  // finally multiply by inverse mass matrix
879  Ulam = invMass*F;
880 
881  // fill column of matrix (Qmat is in column major format)
882  Vmath::Vcopy(ncoeffs,&ulam[0],1,&(Qmat.GetPtr())[0]+j*ncoeffs,1);
883  }
884  }
885  break;
886  // Matrix K (P23)
887  case StdRegions::eHybridDGHelmBndLam:
888  {
889  int i,j,e,cnt;
890  int order_e, nquad_e;
891  int nbndry = NumDGBndryCoeffs();
892  int coordim = GetCoordim();
893  int nedges = GetNedges();
894  NekDouble tau = mkey.GetConstFactor(StdRegions::eFactorTau);
895 
896  Array<OneD,NekDouble> work, varcoeff_work;
897  Array<OneD,const Array<OneD, NekDouble> > normals;
898  Array<OneD,StdRegions::StdExpansionSharedPtr> EdgeExp(nedges);
899  Array<OneD, NekDouble> lam(nbndry);
900 
901  Array<OneD,unsigned int> emap;
902  Array<OneD, int> sign;
903  StdRegions::Orientation edgedir;
904 
905  // declare matrix space
906  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(nbndry, nbndry);
907  DNekMat &BndMat = *returnval;
908 
909  DNekScalMatSharedPtr LamToQ[3];
910 
911  // Matrix to map Lambda to U
912  MatrixKey LamToUkey(StdRegions::eHybridDGLamToU, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
913  DNekScalMat &LamToU = *GetLocMatrix(LamToUkey);
914 
915  // Matrix to map Lambda to Q0
916  MatrixKey LamToQ0key(StdRegions::eHybridDGLamToQ0, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
917  LamToQ[0] = GetLocMatrix(LamToQ0key);
918 
919  // Matrix to map Lambda to Q1
920  MatrixKey LamToQ1key(StdRegions::eHybridDGLamToQ1, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
921  LamToQ[1] = GetLocMatrix(LamToQ1key);
922 
923  // Matrix to map Lambda to Q2 for 3D coordinates
924  if (coordim == 3)
925  {
926  MatrixKey LamToQ2key(StdRegions::eHybridDGLamToQ2, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
927  LamToQ[2] = GetLocMatrix(LamToQ2key);
928  }
929 
930  // Set up edge segment expansions from local geom info
931  for(i = 0; i < nedges; ++i)
932  {
933  EdgeExp[i] = GetEdgeExp(i);
934  }
935 
936  // Set up matrix derived from <mu, Q_lam.n - \tau (U_lam - Lam) >
937  for(i = 0; i < nbndry; ++i)
938  {
939  cnt = 0;
940 
941  Vmath::Zero(nbndry,lam,1);
942  lam[i] = 1.0;
943  SetTraceToGeomOrientation(EdgeExp,lam);
944 
945  for(e = 0; e < nedges; ++e)
946  {
947  order_e = EdgeExp[e]->GetNcoeffs();
948  nquad_e = EdgeExp[e]->GetNumPoints(0);
949 
950  normals = GetEdgeNormal(e);
951  edgedir = GetEorient(e);
952 
953  work = Array<OneD,NekDouble>(nquad_e);
954  varcoeff_work = Array<OneD, NekDouble>(nquad_e);
955 
956  GetEdgeToElementMap(e,edgedir,emap,sign);
957 
958 
959  StdRegions::VarCoeffType VarCoeff[3] = {StdRegions::eVarCoeffD00,
960  StdRegions::eVarCoeffD11,
961  StdRegions::eVarCoeffD22};
962  const StdRegions::VarCoeffMap &varcoeffs = mkey.GetVarCoeffs();
963  StdRegions::VarCoeffMap::const_iterator x;
964 
965  // Q0 * n0 (BQ_0 terms)
966  Array<OneD, NekDouble> edgeCoeffs(order_e);
967  Array<OneD, NekDouble> edgePhys (nquad_e);
968  for(j = 0; j < order_e; ++j)
969  {
970  edgeCoeffs[j] = sign[j]*(*LamToQ[0])(emap[j],i);
971  }
972 
973  EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
974 // @TODO Var coeffs
975  // Multiply by variable coefficient
976 // if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
977 // {
978 // GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
979 // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
980 // }
981 
982  Vmath::Vmul(nquad_e, normals[0], 1, edgePhys, 1, work, 1);
983 
984  // Q1 * n1 (BQ_1 terms)
985  for(j = 0; j < order_e; ++j)
986  {
987  edgeCoeffs[j] = sign[j]*(*LamToQ[1])(emap[j],i);
988  }
989 
990  EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
991 
992 // @TODO var coeffs
993  // Multiply by variable coefficients
994 // if ((x = varcoeffs.find(VarCoeff[1])) != varcoeffs.end())
995 // {
996 // GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
997 // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
998 // }
999 
1000  Vmath::Vvtvp(nquad_e, normals[1], 1, edgePhys, 1,
1001  work, 1, work, 1);
1002 
1003  // Q2 * n2 (BQ_2 terms)
1004  if (coordim == 3)
1005  {
1006  for(j = 0; j < order_e; ++j)
1007  {
1008  edgeCoeffs[j] = sign[j]*(*LamToQ[2])(emap[j],i);
1009  }
1010 
1011  EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
1012 // @TODO var coeffs
1013  // Multiply by variable coefficients
1014 // if ((x = varcoeffs.find(VarCoeff[2])) != varcoeffs.end())
1015 // {
1016 // GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
1017 // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
1018 // }
1019 
1020  Vmath::Vvtvp(nquad_e, normals[2], 1, edgePhys, 1,
1021  work, 1, work, 1);
1022  }
1023 
1024  if (m_negatedNormals[e])
1025  {
1026  Vmath::Neg(nquad_e, work, 1);
1027  }
1028 
1029  // - tau (ulam - lam)
1030  // Corresponds to the G and BU terms.
1031  for(j = 0; j < order_e; ++j)
1032  {
1033  edgeCoeffs[j] = sign[j]*LamToU(emap[j],i) - lam[cnt+j];
1034  }
1035 
1036  EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
1037 
1038  // Multiply by variable coefficients
1039  if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
1040  {
1041  GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
1042  Vmath::Vmul(nquad_e,varcoeff_work,1,edgePhys,1,edgePhys,1);
1043  }
1044 
1045  Vmath::Svtvp(nquad_e,-tau,edgePhys,1,
1046  work,1,work,1);
1047 /// TODO: Add variable coeffs
1048  EdgeExp[e]->IProductWRTBase(work, edgeCoeffs);
1049 
1050  EdgeExp[e]->SetCoeffsToOrientation(edgeCoeffs, edgedir);
1051 
1052  for(j = 0; j < order_e; ++j)
1053  {
1054  BndMat(cnt+j,i) = edgeCoeffs[j];
1055  }
1056 
1057  cnt += order_e;
1058  }
1059  }
1060  }
1061  break;
1062  //HDG postprocessing
1063  case StdRegions::eInvLaplacianWithUnityMean:
1064  {
1065  MatrixKey lapkey(StdRegions::eLaplacian, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1066  DNekScalMat &LapMat = *GetLocMatrix(lapkey);
1067 
1068  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(LapMat.GetRows(),LapMat.GetColumns());
1069  DNekMatSharedPtr lmat = returnval;
1070 
1071  (*lmat) = LapMat;
1072 
1073  // replace first column with inner product wrt 1
1074  int nq = GetTotPoints();
1075  Array<OneD, NekDouble> tmp(nq);
1076  Array<OneD, NekDouble> outarray(m_ncoeffs);
1077  Vmath::Fill(nq,1.0,tmp,1);
1078  IProductWRTBase(tmp, outarray);
1079 
1080  Vmath::Vcopy(m_ncoeffs,&outarray[0],1,
1081  &(lmat->GetPtr())[0],1);
1082 
1083  lmat->Invert();
1084  }
1085  break;
1086  default:
1087  ASSERTL0(false,"This matrix type cannot be generated from this class");
1088  break;
1089  }
1090 
1091  return returnval;
1092  }
1093 
1094  //Evaluate Coefficients of weak deriviative in the direction dir
1095  //given the input coefficicents incoeffs and the imposed
1096  //boundary values in EdgeExp (which will have its phys space updated);
1097  void Expansion2D::v_DGDeriv(
1098  int dir,
1099  const Array<OneD, const NekDouble> &incoeffs,
1100  Array<OneD,StdRegions::StdExpansionSharedPtr> &EdgeExp,
1101  Array<OneD, Array<OneD, NekDouble> > &edgeCoeffs,
1102  Array<OneD, NekDouble> &out_d)
1103  {
1104  StdRegions::MatrixType DerivType[3] = {StdRegions::eWeakDeriv0,
1105  StdRegions::eWeakDeriv1,
1106  StdRegions::eWeakDeriv2};
1107 
1108  int ncoeffs = GetNcoeffs();
1109 
1110  DNekScalMat &InvMass = *GetLocMatrix(StdRegions::eInvMass);
1111  DNekScalMat &Dmat = *GetLocMatrix(DerivType[dir]);
1112 
1113  Array<OneD, NekDouble> coeffs = incoeffs;
1114  DNekVec Coeffs (ncoeffs,coeffs, eWrapper);
1115 
1116  Coeffs = Transpose(Dmat)*Coeffs;
1117  Vmath::Neg(ncoeffs, coeffs,1);
1118 
1119  // Add the boundary integral including the relevant part of
1120  // the normal
1121  AddNormTraceInt(dir, EdgeExp, edgeCoeffs, coeffs);
1122 
1123  DNekVec Out_d (ncoeffs,out_d,eWrapper);
1124 
1125  Out_d = InvMass*Coeffs;
1126  }
1127 
1128  enum BndToLocMatrixMapType
1129  {
1130  eBndToFullMatrixCG,
1131  eBndToBndMatrixCG,
1132  eBndToTraceMatrixDG
1133  };
1134 
1135  void Expansion2D::v_AddRobinMassMatrix(const int edge, const Array<OneD, const NekDouble > &primCoeffs, DNekMatSharedPtr &inoutmat)
1136  {
1137  ASSERTL1(IsBoundaryInteriorExpansion(),
1138  "Not set up for non boundary-interior expansions");
1139  ASSERTL1(inoutmat->GetRows() == inoutmat->GetColumns(),
1140  "Assuming that input matrix was square");
1141  int i,j;
1142  int id1,id2;
1143  ExpansionSharedPtr edgeExp = m_edgeExp[edge].lock();
1144  int order_e = edgeExp->GetNcoeffs();
1145 
1146  Array<OneD,unsigned int> map;
1147  Array<OneD,int> sign;
1148 
1149  StdRegions::VarCoeffMap varcoeffs;
1150  varcoeffs[StdRegions::eVarCoeffMass] = primCoeffs;
1151 
1152  LocalRegions::MatrixKey mkey(StdRegions::eMass,LibUtilities::eSegment, *edgeExp, StdRegions::NullConstFactorMap, varcoeffs);
1153  DNekScalMat &edgemat = *edgeExp->GetLocMatrix(mkey);
1154 
1155  // Now need to identify a map which takes the local edge
1156  // mass matrix to the matrix stored in inoutmat;
1157  // This can currently be deduced from the size of the matrix
1158 
1159  // - if inoutmat.m_rows() == v_NCoeffs() it is a full
1160  // matrix system
1161 
1162  // - if inoutmat.m_rows() == v_NumBndCoeffs() it is a
1163  // boundary CG system
1164 
1165  // - if inoutmat.m_rows() == v_NumDGBndCoeffs() it is a
1166  // trace DG system
1167  int rows = inoutmat->GetRows();
1168 
1169  if (rows == GetNcoeffs())
1170  {
1171  GetEdgeToElementMap(edge,v_GetEorient(edge),map,sign);
1172  }
1173  else if(rows == NumBndryCoeffs())
1174  {
1175  int nbndry = NumBndryCoeffs();
1176  Array<OneD,unsigned int> bmap(nbndry);
1177 
1178  GetEdgeToElementMap(edge,v_GetEorient(edge),map,sign);
1179 
1180  GetBoundaryMap(bmap);
1181 
1182  for(i = 0; i < order_e; ++i)
1183  {
1184  for(j = 0; j < nbndry; ++j)
1185  {
1186  if(map[i] == bmap[j])
1187  {
1188  map[i] = j;
1189  break;
1190  }
1191  }
1192  ASSERTL1(j != nbndry,"Did not find number in map");
1193  }
1194  }
1195  else if (rows == NumDGBndryCoeffs())
1196  {
1197  // possibly this should be a separate method
1198  int cnt = 0;
1199  map = Array<OneD, unsigned int> (order_e);
1200  sign = Array<OneD, int> (order_e,1);
1201 
1202  for(i = 0; i < edge; ++i)
1203  {
1204  cnt += GetEdgeNcoeffs(i);
1205  }
1206 
1207  for(i = 0; i < order_e; ++i)
1208  {
1209  map[i] = cnt++;
1210  }
1211  // check for mapping reversal
1212  if(GetEorient(edge) == StdRegions::eBackwards)
1213  {
1214  switch(edgeExp->GetBasis(0)->GetBasisType())
1215  {
1216  case LibUtilities::eGauss_Lagrange:
1217  reverse( map.get() , map.get()+order_e);
1218  break;
1219  case LibUtilities::eGLL_Lagrange:
1220  reverse( map.get() , map.get()+order_e);
1221  break;
1222  case LibUtilities::eModified_A:
1223  {
1224  swap(map[0],map[1]);
1225  for(i = 3; i < order_e; i+=2)
1226  {
1227  sign[i] = -1;
1228  }
1229  }
1230  break;
1231  default:
1232  ASSERTL0(false,"Edge boundary type not valid for this method");
1233  }
1234  }
1235  }
1236  else
1237  {
1238  ASSERTL0(false,"Could not identify matrix type from dimension");
1239  }
1240 
1241  for(i = 0; i < order_e; ++i)
1242  {
1243  id1 = map[i];
1244  for(j = 0; j < order_e; ++j)
1245  {
1246  id2 = map[j];
1247  (*inoutmat)(id1,id2) += edgemat(i,j)*sign[i]*sign[j];
1248  }
1249  }
1250  }
1251 
1252  /**
1253  * Given an edge and vector of element coefficients:
1254  * - maps those elemental coefficients corresponding to the edge into
1255  * an edge-vector.
1256  * - resets the element coefficients
1257  * - multiplies the edge vector by the edge mass matrix
1258  * - maps the edge coefficients back onto the elemental coefficients
1259  */
1260  void Expansion2D::v_AddRobinEdgeContribution(const int edgeid, const Array<OneD, const NekDouble> &primCoeffs, Array<OneD, NekDouble> &coeffs)
1261  {
1262  ASSERTL1(IsBoundaryInteriorExpansion(),
1263  "Not set up for non boundary-interior expansions");
1264  int i;
1265  ExpansionSharedPtr edgeExp = m_edgeExp[edgeid].lock();
1266  int order_e = edgeExp->GetNcoeffs();
1267 
1268  Array<OneD,unsigned int> map;
1269  Array<OneD,int> sign;
1270 
1271  StdRegions::VarCoeffMap varcoeffs;
1272  varcoeffs[StdRegions::eVarCoeffMass] = primCoeffs;
1273 
1274  LocalRegions::MatrixKey mkey(StdRegions::eMass,LibUtilities::eSegment, *edgeExp, StdRegions::NullConstFactorMap, varcoeffs);
1275  DNekScalMat &edgemat = *edgeExp->GetLocMatrix(mkey);
1276 
1277  NekVector<NekDouble> vEdgeCoeffs (order_e);
1278 
1279  GetEdgeToElementMap(edgeid,v_GetEorient(edgeid),map,sign);
1280 
1281  for (i = 0; i < order_e; ++i)
1282  {
1283  vEdgeCoeffs[i] = coeffs[map[i]]*sign[i];
1284  }
1285  Vmath::Zero(GetNcoeffs(), coeffs, 1);
1286 
1287  vEdgeCoeffs = edgemat * vEdgeCoeffs;
1288 
1289  for (i = 0; i < order_e; ++i)
1290  {
1291  coeffs[map[i]] = vEdgeCoeffs[i]*sign[i];
1292  }
1293  }
1294 
1295  DNekMatSharedPtr Expansion2D::v_BuildVertexMatrix(
1296  const DNekScalMatSharedPtr &r_bnd)
1297  {
1298  MatrixStorage storage = eFULL;
1299  DNekMatSharedPtr m_vertexmatrix;
1300 
1301  int nVerts, vid1, vid2, vMap1, vMap2;
1302  NekDouble VertexValue;
1303 
1304  nVerts = GetNverts();
1305 
1306  m_vertexmatrix =
1307  MemoryManager<DNekMat>::AllocateSharedPtr(
1308  nVerts, nVerts, 0.0, storage);
1309  DNekMat &VertexMat = (*m_vertexmatrix);
1310 
1311  for (vid1 = 0; vid1 < nVerts; ++vid1)
1312  {
1313  vMap1 = GetVertexMap(vid1);
1314 
1315  for (vid2 = 0; vid2 < nVerts; ++vid2)
1316  {
1317  vMap2 = GetVertexMap(vid2);
1318  VertexValue = (*r_bnd)(vMap1, vMap2);
1319  VertexMat.SetValue(vid1, vid2, VertexValue);
1320  }
1321  }
1322 
1323  return m_vertexmatrix;
1324  }
1325 
1326  Array<OneD, unsigned int> Expansion2D::v_GetEdgeInverseBoundaryMap(
1327  int eid)
1328  {
1329  int n, j;
1330  int nEdgeCoeffs;
1331  int nBndCoeffs = NumBndryCoeffs();
1332 
1333  Array<OneD, unsigned int> bmap(nBndCoeffs);
1334  GetBoundaryMap(bmap);
1335 
1336  // Map from full system to statically condensed system (i.e reverse
1337  // GetBoundaryMap)
1338  map<int, int> invmap;
1339  for (j = 0; j < nBndCoeffs; ++j)
1340  {
1341  invmap[bmap[j]] = j;
1342  }
1343 
1344  // Number of interior edge coefficients
1345  nEdgeCoeffs = GetEdgeNcoeffs(eid) - 2;
1346 
1347  const SpatialDomains::Geometry2DSharedPtr &geom = GetGeom2D();
1348 
1349  Array<OneD, unsigned int> edgemaparray(nEdgeCoeffs);
1350  Array<OneD, unsigned int> maparray (nEdgeCoeffs);
1351  Array<OneD, int> signarray (nEdgeCoeffs, 1);
1352  StdRegions::Orientation eOrient = geom->GetEorient(eid);
1353 
1354  // maparray is the location of the edge within the matrix
1355  GetEdgeInteriorMap(eid, eOrient, maparray, signarray);
1356 
1357  for (n = 0; n < nEdgeCoeffs; ++n)
1358  {
1359  edgemaparray[n] = invmap[maparray[n]];
1360  }
1361 
1362  return edgemaparray;
1363  }
1364 
1365  void Expansion2D::v_SetUpPhysNormals(const int edge)
1366  {
1367  ComputeEdgeNormal(edge);
1368  }
1369 
1370  const StdRegions::NormalVector &Expansion2D::v_GetEdgeNormal(const int edge) const
1371  {
1372  std::map<int, StdRegions::NormalVector>::const_iterator x;
1373  x = m_edgeNormals.find(edge);
1374  ASSERTL0 (x != m_edgeNormals.end(),
1375  "Edge normal not computed.");
1376  return x->second;
1377  }
1378 
1379  const StdRegions::NormalVector &Expansion2D::v_GetSurfaceNormal(
1380  const int id) const
1381  {
1382  return v_GetEdgeNormal(id);
1383  }
1384 
1385  void Expansion2D::v_NegateEdgeNormal(const int edge)
1386  {
1387  m_negatedNormals[edge] = true;
1388  for (int i = 0; i < GetCoordim(); ++i)
1389  {
1390  Vmath::Neg(m_edgeNormals[edge][i].num_elements(),
1391  m_edgeNormals[edge][i], 1);
1392  }
1393  }
1394 
1395  bool Expansion2D::v_EdgeNormalNegated(const int edge)
1396  {
1397  return m_negatedNormals[edge];
1398  }
1399  }
1400 }
Generated on Sun Mar 15 2015 00:11:31 for Nektar++ by   doxygen 1.8.1.2