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