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