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