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