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