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