Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A,
1107  "Method set up only for modified modal bases curretly");
1108 
1109  int i,j;
1110  int id1,id2;
1111  Expansion2DSharedPtr faceExp = m_faceExp[face].lock();
1112  int order_f = faceExp->GetNcoeffs();
1113 
1114  Array<OneD, unsigned int> map;
1115  Array<OneD, int> sign;
1116 
1117  StdRegions::VarCoeffMap varcoeffs;
1118  varcoeffs[StdRegions::eVarCoeffMass] = primCoeffs;
1119 
1120  LibUtilities::ShapeType shapeType =
1121  faceExp->DetShapeType();
1122 
1123  LocalRegions::MatrixKey mkey(StdRegions::eMass,
1124  shapeType,
1125  *faceExp,
1126  StdRegions::NullConstFactorMap,
1127  varcoeffs);
1128 
1129  DNekScalMat &facemat = *faceExp->GetLocMatrix(mkey);
1130 
1131  // Now need to identify a map which takes the local face
1132  // mass matrix to the matrix stored in inoutmat;
1133  // This can currently be deduced from the size of the matrix
1134 
1135  // - if inoutmat.m_rows() == v_NCoeffs() it is a full
1136  // matrix system
1137 
1138  // - if inoutmat.m_rows() == v_GetNverts() it is a vertex space
1139  // preconditioner.
1140 
1141  // - if inoutmat.m_rows() == v_NumBndCoeffs() it is a
1142  // boundary CG system
1143 
1144  // - if inoutmat.m_rows() == v_NumDGBndCoeffs() it is a
1145  // trace DG system; still needs implementing.
1146  int rows = inoutmat->GetRows();
1147 
1148  if (rows == GetNcoeffs())
1149  {
1150  GetFaceToElementMap(face,GetForient(face),map,sign);
1151  }
1152  else if (rows == GetNverts())
1153  {
1154  int nfvert = faceExp->GetNverts();
1155 
1156  // Need to find where linear vertices are in facemat
1157  Array<OneD, unsigned int> linmap;
1158  Array<OneD, int> linsign;
1159 
1160  // Use a linear expansion to get correct mapping
1161  GetLinStdExp()->GetFaceToElementMap(face,GetForient(face),linmap, linsign);
1162 
1163  // zero out sign map to remove all other modes
1164  sign = Array<OneD, int> (order_f,0);
1165  map = Array<OneD, unsigned int>(order_f,(unsigned int)0);
1166 
1167  int fmap;
1168  // Reset sign map to only have contribution from vertices
1169  for(i = 0; i < nfvert; ++i)
1170  {
1171  fmap = faceExp->GetVertexMap(i,true);
1172  sign[fmap] = 1;
1173 
1174  // need to reset map
1175  map[fmap] = linmap[i];
1176  }
1177  }
1178  else if(rows == NumBndryCoeffs())
1179  {
1180  int nbndry = NumBndryCoeffs();
1181  Array<OneD,unsigned int> bmap(nbndry);
1182 
1183  GetFaceToElementMap(face,GetForient(face),map,sign);
1184  GetBoundaryMap(bmap);
1185 
1186  for(i = 0; i < order_f; ++i)
1187  {
1188  for(j = 0; j < nbndry; ++j)
1189  {
1190  if(map[i] == bmap[j])
1191  {
1192  map[i] = j;
1193  break;
1194  }
1195  }
1196  ASSERTL1(j != nbndry,"Did not find number in map");
1197  }
1198  }
1199  else if (rows == NumDGBndryCoeffs())
1200  {
1201  // possibly this should be a separate method
1202  int cnt = 0;
1203  map = Array<OneD, unsigned int> (order_f);
1204  sign = Array<OneD, int> (order_f,1);
1205 
1206  StdRegions::IndexMapKey ikey1(
1207  StdRegions::eFaceToElement, DetShapeType(),
1208  GetBasisNumModes(0), GetBasisNumModes(1), GetBasisNumModes(2),
1209  face, GetForient(face));
1210  StdRegions::IndexMapValuesSharedPtr map1 =
1211  StdExpansion::GetIndexMap(ikey1);
1212  StdRegions::IndexMapKey ikey2(
1213  StdRegions::eFaceToElement,
1214  DetShapeType(),
1215  GetBasisNumModes(0),
1216  GetBasisNumModes(1),
1217  GetBasisNumModes(2),
1218  face,
1219  StdRegions::eDir1FwdDir1_Dir2FwdDir2);
1220  StdRegions::IndexMapValuesSharedPtr map2 =
1221  StdExpansion::GetIndexMap(ikey2);
1222 
1223  ASSERTL1((*map1).num_elements() == (*map2).num_elements(),
1224  "There is an error with the GetFaceToElementMap");
1225 
1226  for (i = 0; i < face; ++i)
1227  {
1228  cnt += GetFaceNcoeffs(i);
1229  }
1230 
1231  for(i = 0; i < (*map1).num_elements(); ++i)
1232  {
1233  int idx = -1;
1234 
1235  for(j = 0; j < (*map2).num_elements(); ++j)
1236  {
1237  if((*map1)[i].index == (*map2)[j].index)
1238  {
1239  idx = j;
1240  break;
1241  }
1242  }
1243 
1244  ASSERTL2(idx >= 0, "Index not found");
1245  map [i] = idx + cnt;
1246  sign[i] = (*map2)[idx].sign;
1247  }
1248  }
1249  else
1250  {
1251  ASSERTL0(false,"Could not identify matrix type from dimension");
1252  }
1253 
1254  for(i = 0; i < order_f; ++i)
1255  {
1256  id1 = map[i];
1257  for(j = 0; j < order_f; ++j)
1258  {
1259  id2 = map[j];
1260  (*inoutmat)(id1,id2) += facemat(i,j)*sign[i]*sign[j];
1261  }
1262  }
1263  }
1264 
1265  DNekMatSharedPtr Expansion3D::v_BuildVertexMatrix(
1266  const DNekScalMatSharedPtr &r_bnd)
1267  {
1268  MatrixStorage storage = eFULL;
1269  DNekMatSharedPtr m_vertexmatrix;
1270 
1271  int nVerts, vid1, vid2, vMap1, vMap2;
1272  NekDouble VertexValue;
1273 
1274  nVerts = GetNverts();
1275 
1276  m_vertexmatrix =
1277  MemoryManager<DNekMat>::AllocateSharedPtr(
1278  nVerts, nVerts, 0.0, storage);
1279  DNekMat &VertexMat = (*m_vertexmatrix);
1280 
1281  for (vid1 = 0; vid1 < nVerts; ++vid1)
1282  {
1283  vMap1 = GetVertexMap(vid1,true);
1284 
1285  for (vid2 = 0; vid2 < nVerts; ++vid2)
1286  {
1287  vMap2 = GetVertexMap(vid2,true);
1288  VertexValue = (*r_bnd)(vMap1, vMap2);
1289  VertexMat.SetValue(vid1, vid2, VertexValue);
1290  }
1291  }
1292 
1293  return m_vertexmatrix;
1294  }
1295 
1296  DNekMatSharedPtr Expansion3D::v_BuildTransformationMatrix(
1297  const DNekScalMatSharedPtr &r_bnd,
1298  const StdRegions::MatrixType matrixType)
1299  {
1300  int nVerts, nEdges;
1301  int eid, fid, vid, n, i;
1302 
1303  int nBndCoeffs = NumBndryCoeffs();
1304 
1305  const SpatialDomains::Geometry3DSharedPtr &geom = GetGeom3D();
1306 
1307  // Get geometric information about this element
1308  nVerts = GetNverts();
1309  nEdges = GetNedges();
1310 
1311  /*************************************/
1312  /* Vetex-edge & vertex-face matrices */
1313  /*************************************/
1314 
1315  /**
1316  * The matrix component of \f$\mathbf{R}\f$ is given by \f[
1317  * \mathbf{R^{T}_{v}}=
1318  * -\mathbf{S}^{-1}_{ef,ef}\mathbf{S}^{T}_{v,ef}\f]
1319  *
1320  * For every vertex mode we extract the submatrices from statically
1321  * condensed matrix \f$\mathbf{S}\f$ corresponding to the coupling
1322  * between the attached edges and faces of a vertex
1323  * (\f$\mathbf{S_{ef,ef}}\f$). This matrix is then inverted and
1324  * multiplied by the submatrix representing the coupling between a
1325  * vertex and the attached edges and faces
1326  * (\f$\mathbf{S_{v,ef}}\f$).
1327  */
1328 
1329  int nmodes;
1330  int m;
1331  NekDouble VertexEdgeFaceValue;
1332 
1333  // The number of connected edges/faces is 3 (for all elements)
1334  int nConnectedEdges = 3;
1335  int nConnectedFaces = 3;
1336 
1337  // Location in the matrix
1338  Array<OneD, Array<OneD, unsigned int> >
1339  MatEdgeLocation(nConnectedEdges);
1340  Array<OneD, Array<OneD, unsigned int> >
1341  MatFaceLocation(nConnectedFaces);
1342 
1343  // Define storage for vertex transpose matrix and zero all entries
1344  MatrixStorage storage = eFULL;
1345  DNekMatSharedPtr m_transformationmatrix;
1346  DNekMatSharedPtr m_transposedtransformationmatrix;
1347 
1348  m_transformationmatrix =
1349  MemoryManager<DNekMat>::AllocateSharedPtr(
1350  nBndCoeffs, nBndCoeffs, 0.0, storage);
1351  m_transposedtransformationmatrix =
1352  MemoryManager<DNekMat>::AllocateSharedPtr(
1353  nBndCoeffs, nBndCoeffs, 0.0, storage);
1354 
1355  DNekMat &R = (*m_transformationmatrix);
1356  DNekMat &RT = (*m_transposedtransformationmatrix);
1357 
1358  // Build the vertex-edge/face transform matrix: This matrix is
1359  // constructed from the submatrices corresponding to the couping
1360  // between each vertex and the attached edges/faces
1361  for (vid = 0; vid < nVerts; ++vid)
1362  {
1363  // Row and column size of the vertex-edge/face matrix
1364  int efRow =
1365  GetEdgeNcoeffs (geom->GetVertexEdgeMap(vid, 0)) +
1366  GetEdgeNcoeffs (geom->GetVertexEdgeMap(vid, 1)) +
1367  GetEdgeNcoeffs (geom->GetVertexEdgeMap(vid, 2)) +
1368  GetFaceIntNcoeffs(geom->GetVertexFaceMap(vid, 0)) +
1369  GetFaceIntNcoeffs(geom->GetVertexFaceMap(vid, 1)) +
1370  GetFaceIntNcoeffs(geom->GetVertexFaceMap(vid, 2)) - 6;
1371 
1372  int nedgemodesconnected =
1373  GetEdgeNcoeffs (geom->GetVertexEdgeMap(vid, 0)) +
1374  GetEdgeNcoeffs (geom->GetVertexEdgeMap(vid, 1)) +
1375  GetEdgeNcoeffs (geom->GetVertexEdgeMap(vid, 2)) - 6;
1376  Array<OneD, unsigned int> edgemodearray(nedgemodesconnected);
1377 
1378  int nfacemodesconnected =
1379  GetFaceIntNcoeffs(geom->GetVertexFaceMap(vid, 0)) +
1380  GetFaceIntNcoeffs(geom->GetVertexFaceMap(vid, 1)) +
1381  GetFaceIntNcoeffs(geom->GetVertexFaceMap(vid, 2));
1382  Array<OneD, unsigned int> facemodearray(nfacemodesconnected);
1383 
1384  int offset = 0;
1385 
1386  // Create array of edge modes
1387  for (eid = 0; eid < nConnectedEdges; ++eid)
1388  {
1389  MatEdgeLocation[eid] = GetEdgeInverseBoundaryMap(
1390  geom->GetVertexEdgeMap(vid, eid));
1391  nmodes = MatEdgeLocation[eid].num_elements();
1392 
1393  if (nmodes)
1394  {
1395  Vmath::Vcopy(nmodes, &MatEdgeLocation[eid][0],
1396  1, &edgemodearray[offset], 1);
1397  }
1398 
1399  offset += nmodes;
1400  }
1401 
1402  offset = 0;
1403 
1404  // Create array of face modes
1405  for (fid = 0; fid < nConnectedFaces; ++fid)
1406  {
1407  MatFaceLocation[fid] = GetFaceInverseBoundaryMap(
1408  geom->GetVertexFaceMap(vid, fid));
1409  nmodes = MatFaceLocation[fid].num_elements();
1410 
1411  if (nmodes)
1412  {
1413  Vmath::Vcopy(nmodes, &MatFaceLocation[fid][0],
1414  1, &facemodearray[offset], 1);
1415  }
1416  offset += nmodes;
1417  }
1418 
1419  DNekMatSharedPtr m_vertexedgefacetransformmatrix =
1420  MemoryManager<DNekMat>::AllocateSharedPtr(
1421  1, efRow, 0.0, storage);
1422  DNekMat &Sveft = (*m_vertexedgefacetransformmatrix);
1423 
1424  DNekMatSharedPtr m_vertexedgefacecoupling =
1425  MemoryManager<DNekMat>::AllocateSharedPtr(
1426  1, efRow, 0.0, storage);
1427  DNekMat &Svef = (*m_vertexedgefacecoupling);
1428 
1429  // Vertex-edge coupling
1430  for (n = 0; n < nedgemodesconnected; ++n)
1431  {
1432  // Matrix value for each coefficient location
1433  VertexEdgeFaceValue = (*r_bnd)(GetVertexMap(vid),
1434  edgemodearray[n]);
1435 
1436  // Set the value in the vertex edge/face matrix
1437  Svef.SetValue(0, n, VertexEdgeFaceValue);
1438  }
1439 
1440  // Vertex-face coupling
1441  for (n = 0; n < nfacemodesconnected; ++n)
1442  {
1443  // Matrix value for each coefficient location
1444  VertexEdgeFaceValue = (*r_bnd)(GetVertexMap(vid),
1445  facemodearray[n]);
1446 
1447  // Set the value in the vertex edge/face matrix
1448  Svef.SetValue(0, n + nedgemodesconnected, VertexEdgeFaceValue);
1449  }
1450 
1451  /*
1452  * Build the edge-face transform matrix: This matrix is
1453  * constructed from the submatrices corresponding to the couping
1454  * between the edges and faces on the attached faces/edges of a
1455  * vertex.
1456  */
1457 
1458  // Allocation of matrix to store edge/face-edge/face coupling
1459  DNekMatSharedPtr m_edgefacecoupling =
1460  MemoryManager<DNekMat>::AllocateSharedPtr(
1461  efRow, efRow, 0.0, storage);
1462  DNekMat &Sefef = (*m_edgefacecoupling);
1463 
1464  NekDouble EdgeEdgeValue, FaceFaceValue;
1465 
1466  // Edge-edge coupling (S_{ee})
1467  for (m = 0; m < nedgemodesconnected; ++m)
1468  {
1469  for (n = 0; n < nedgemodesconnected; ++n)
1470  {
1471  // Matrix value for each coefficient location
1472  EdgeEdgeValue = (*r_bnd)(edgemodearray[n],
1473  edgemodearray[m]);
1474 
1475  // Set the value in the vertex edge/face matrix
1476  Sefef.SetValue(n, m, EdgeEdgeValue);
1477  }
1478  }
1479 
1480  // Face-face coupling (S_{ff})
1481  for (n = 0; n < nfacemodesconnected; ++n)
1482  {
1483  for (m = 0; m < nfacemodesconnected; ++m)
1484  {
1485  // Matrix value for each coefficient location
1486  FaceFaceValue = (*r_bnd)(facemodearray[n],
1487  facemodearray[m]);
1488  // Set the value in the vertex edge/face matrix
1489  Sefef.SetValue(nedgemodesconnected + n,
1490  nedgemodesconnected + m, FaceFaceValue);
1491  }
1492  }
1493 
1494  // Edge-face coupling (S_{ef} and trans(S_{ef}))
1495  for (n = 0; n < nedgemodesconnected; ++n)
1496  {
1497  for (m = 0; m < nfacemodesconnected; ++m)
1498  {
1499  // Matrix value for each coefficient location
1500  FaceFaceValue = (*r_bnd)(edgemodearray[n],
1501  facemodearray[m]);
1502 
1503  // Set the value in the vertex edge/face matrix
1504  Sefef.SetValue(n,
1505  nedgemodesconnected + m,
1506  FaceFaceValue);
1507 
1508  // and transpose
1509  Sefef.SetValue(nedgemodesconnected + m,
1510  n,
1511  FaceFaceValue);
1512  }
1513  }
1514 
1515 
1516  // Invert edge-face coupling matrix
1517  if (efRow)
1518  {
1519  Sefef.Invert();
1520 
1521  //R_{v}=-S_{v,ef}inv(S_{ef,ef})
1522  Sveft = -Svef * Sefef;
1523  }
1524 
1525  // Populate R with R_{ve} components
1526  for (n = 0; n < edgemodearray.num_elements(); ++n)
1527  {
1528  RT.SetValue(edgemodearray[n], GetVertexMap(vid),
1529  Sveft(0, n));
1530  R.SetValue(GetVertexMap(vid), edgemodearray[n],
1531  Sveft(0, n));
1532  }
1533 
1534  // Populate R with R_{vf} components
1535  for (n = 0; n < facemodearray.num_elements(); ++n)
1536  {
1537  RT.SetValue(facemodearray[n], GetVertexMap(vid),
1538  Sveft(0, n + nedgemodesconnected));
1539  R.SetValue(GetVertexMap(vid), facemodearray[n],
1540  Sveft(0, n + nedgemodesconnected));
1541  }
1542  }
1543 
1544  /********************/
1545  /* edge-face matrix */
1546  /********************/
1547 
1548  /*
1549  * The matrix component of \f$\mathbf{R}\f$ is given by \f[
1550  * \mathbf{R^{T}_{ef}}=-\mathbf{S}^{-1}_{ff}\mathbf{S}^{T}_{ef}\f]
1551  *
1552  * For each edge extract the submatrices from statically condensed
1553  * matrix \f$\mathbf{S}\f$ corresponding to inner products of modes
1554  * on the two attached faces within themselves as well as the
1555  * coupling matrix between the two faces
1556  * (\f$\mathbf{S}_{ff}\f$). This matrix of face coupling is then
1557  * inverted and multiplied by the submatrices of corresponding to
1558  * the coupling between the edge and attached faces
1559  * (\f$\mathbf{S}_{ef}\f$).
1560  */
1561 
1562  NekDouble EdgeFaceValue, FaceFaceValue;
1563  int efCol, efRow, nedgemodes;
1564 
1565  // Number of attached faces is always 2
1566  nConnectedFaces = 2;
1567 
1568  // Location in the matrix
1569  MatEdgeLocation = Array<OneD, Array<OneD, unsigned int> >
1570  (nEdges);
1571  MatFaceLocation = Array<OneD, Array<OneD, unsigned int> >
1572  (nConnectedFaces);
1573 
1574  // Build the edge/face transform matrix: This matrix is constructed
1575  // from the submatrices corresponding to the couping between a
1576  // specific edge and the two attached faces.
1577  for (eid = 0; eid < nEdges; ++eid)
1578  {
1579  // Row and column size of the vertex-edge/face matrix
1580  efCol = GetFaceIntNcoeffs(geom->GetEdgeFaceMap(eid, 0)) +
1581  GetFaceIntNcoeffs(geom->GetEdgeFaceMap(eid, 1));
1582  efRow = GetEdgeNcoeffs(eid) - 2;
1583 
1584  // Edge-face coupling matrix
1585  DNekMatSharedPtr m_efedgefacecoupling =
1586  MemoryManager<DNekMat>::AllocateSharedPtr(
1587  efRow, efCol, 0.0, storage);
1588  DNekMat &Mef = (*m_efedgefacecoupling);
1589 
1590  // Face-face coupling matrix
1591  DNekMatSharedPtr m_effacefacecoupling =
1592  MemoryManager<DNekMat>::AllocateSharedPtr(
1593  efCol, efCol, 0.0, storage);
1594  DNekMat &Meff = (*m_effacefacecoupling);
1595 
1596  // Edge-face transformation matrix
1597  DNekMatSharedPtr m_edgefacetransformmatrix =
1598  MemoryManager<DNekMat>::AllocateSharedPtr(
1599  efRow, efCol, 0.0, storage);
1600  DNekMat &Meft = (*m_edgefacetransformmatrix);
1601 
1602  int nfacemodesconnected =
1603  GetFaceIntNcoeffs(geom->GetEdgeFaceMap(eid, 0)) +
1604  GetFaceIntNcoeffs(geom->GetEdgeFaceMap(eid, 1));
1605  Array<OneD, unsigned int>
1606  facemodearray(nfacemodesconnected);
1607 
1608  // Create array of edge modes
1609  Array<OneD, unsigned int> inedgearray
1610  = GetEdgeInverseBoundaryMap(eid);
1611  nedgemodes = GetEdgeNcoeffs(eid) - 2;
1612  Array<OneD, unsigned int> edgemodearray(nedgemodes);
1613 
1614  if (nedgemodes)
1615  {
1616  Vmath::Vcopy(nedgemodes, &inedgearray[0],
1617  1, &edgemodearray[0], 1);
1618  }
1619 
1620  int offset = 0;
1621 
1622  // Create array of face modes
1623  for (fid = 0; fid < nConnectedFaces; ++fid)
1624  {
1625  MatFaceLocation[fid] = GetFaceInverseBoundaryMap(
1626  geom->GetEdgeFaceMap(eid, fid));
1627  nmodes = MatFaceLocation[fid].num_elements();
1628 
1629  if (nmodes)
1630  {
1631  Vmath::Vcopy(nmodes, &MatFaceLocation[fid][0],
1632  1, &facemodearray[offset], 1);
1633  }
1634  offset += nmodes;
1635  }
1636 
1637  // Edge-face coupling
1638  for (n = 0; n < nedgemodes; ++n)
1639  {
1640  for (m = 0; m < nfacemodesconnected; ++m)
1641  {
1642  // Matrix value for each coefficient location
1643  EdgeFaceValue = (*r_bnd)(edgemodearray[n],
1644  facemodearray[m]);
1645 
1646  // Set the value in the edge/face matrix
1647  Mef.SetValue(n, m, EdgeFaceValue);
1648  }
1649  }
1650 
1651  // Face-face coupling
1652  for (n = 0; n < nfacemodesconnected; ++n)
1653  {
1654  for (m = 0; m < nfacemodesconnected; ++m)
1655  {
1656  // Matrix value for each coefficient location
1657  FaceFaceValue = (*r_bnd)(facemodearray[n],
1658  facemodearray[m]);
1659 
1660  // Set the value in the vertex edge/face matrix
1661  Meff.SetValue(n, m, FaceFaceValue);
1662  }
1663  }
1664 
1665  if (efCol)
1666  {
1667  // Invert edge-face coupling matrix
1668  Meff.Invert();
1669 
1670  // trans(R_{ef})=-S_{ef}*(inv(S_{ff})
1671  Meft = -Mef * Meff;
1672  }
1673 
1674  //Populate transformation matrix with Meft
1675  for (n = 0; n < Meft.GetRows(); ++n)
1676  {
1677  for (m = 0; m < Meft.GetColumns(); ++m)
1678  {
1679  R.SetValue(edgemodearray[n], facemodearray[m],
1680  Meft(n, m));
1681  RT.SetValue(facemodearray[m], edgemodearray[n],
1682  Meft(n, m));
1683  }
1684  }
1685  }
1686 
1687  for (i = 0; i < R.GetRows(); ++i)
1688  {
1689  R.SetValue(i, i, 1.0);
1690  RT.SetValue(i, i, 1.0);
1691  }
1692 
1693  if ((matrixType == StdRegions::ePreconR)||
1694  (matrixType == StdRegions::ePreconRMass))
1695  {
1696  return m_transformationmatrix;
1697  }
1698  else if ((matrixType == StdRegions::ePreconRT)||
1699  (matrixType == StdRegions::ePreconRTMass))
1700  {
1701  return m_transposedtransformationmatrix;
1702  }
1703  else
1704  {
1705  NEKERROR(ErrorUtil::efatal, "unkown matrix type" );
1706  return NullDNekMatSharedPtr;
1707  }
1708  }
1709 
1710  /**
1711  * \brief Build inverse and inverse transposed transformation matrix:
1712  * \f$\mathbf{R^{-1}}\f$ and \f$\mathbf{R^{-T}}\f$
1713  *
1714  * \f\mathbf{R^{-T}}=[\left[\begin{array}{ccc} \mathbf{I} &
1715  * -\mathbf{R}_{ef} & -\mathbf{R}_{ve}+\mathbf{R}_{ve}\mathbf{R}_{vf} \\
1716  * 0 & \mathbf{I} & \mathbf{R}_{ef} \\
1717  * 0 & 0 & \mathbf{I}} \end{array}\right]\f]
1718  */
1719  DNekMatSharedPtr Expansion3D::v_BuildInverseTransformationMatrix(
1720  const DNekScalMatSharedPtr & m_transformationmatrix)
1721  {
1722  int i, j, n, eid = 0, fid = 0;
1723  int nCoeffs = NumBndryCoeffs();
1724  NekDouble MatrixValue;
1725  DNekScalMat &R = (*m_transformationmatrix);
1726 
1727  // Define storage for vertex transpose matrix and zero all entries
1728  MatrixStorage storage = eFULL;
1729 
1730  // Inverse transformation matrix
1731  DNekMatSharedPtr m_inversetransformationmatrix =
1732  MemoryManager<DNekMat>::AllocateSharedPtr(
1733  nCoeffs, nCoeffs, 0.0, storage);
1734  DNekMat &InvR = (*m_inversetransformationmatrix);
1735 
1736  int nVerts = GetNverts();
1737  int nEdges = GetNedges();
1738  int nFaces = GetNfaces();
1739 
1740  int nedgemodes = 0;
1741  int nfacemodes = 0;
1742  int nedgemodestotal = 0;
1743  int nfacemodestotal = 0;
1744 
1745  for (eid = 0; eid < nEdges; ++eid)
1746  {
1747  nedgemodes = GetEdgeNcoeffs(eid) - 2;
1748  nedgemodestotal += nedgemodes;
1749  }
1750 
1751  for (fid = 0; fid < nFaces; ++fid)
1752  {
1753  nfacemodes = GetFaceIntNcoeffs(fid);
1754  nfacemodestotal += nfacemodes;
1755  }
1756 
1757  Array<OneD, unsigned int>
1758  edgemodearray(nedgemodestotal);
1759  Array<OneD, unsigned int>
1760  facemodearray(nfacemodestotal);
1761 
1762  int offset = 0;
1763 
1764  // Create array of edge modes
1765  for (eid = 0; eid < nEdges; ++eid)
1766  {
1767  Array<OneD, unsigned int> edgearray
1768  = GetEdgeInverseBoundaryMap(eid);
1769  nedgemodes = GetEdgeNcoeffs(eid) - 2;
1770 
1771  // Only copy if there are edge modes
1772  if (nedgemodes)
1773  {
1774  Vmath::Vcopy(nedgemodes, &edgearray[0], 1,
1775  &edgemodearray[offset], 1);
1776  }
1777 
1778  offset += nedgemodes;
1779  }
1780 
1781  offset = 0;
1782 
1783  // Create array of face modes
1784  for (fid = 0; fid < nFaces; ++fid)
1785  {
1786  Array<OneD, unsigned int> facearray
1787  = GetFaceInverseBoundaryMap(fid);
1788  nfacemodes = GetFaceIntNcoeffs(fid);
1789 
1790  // Only copy if there are face modes
1791  if (nfacemodes)
1792  {
1793  Vmath::Vcopy(nfacemodes, &facearray[0], 1,
1794  &facemodearray[offset], 1);
1795  }
1796 
1797  offset += nfacemodes;
1798  }
1799 
1800  // Vertex-edge/face
1801  for (i = 0; i < nVerts; ++i)
1802  {
1803  for (j = 0; j < nedgemodestotal; ++j)
1804  {
1805  InvR.SetValue(
1806  GetVertexMap(i), edgemodearray[j],
1807  -R(GetVertexMap(i), edgemodearray[j]));
1808  }
1809  for (j = 0; j < nfacemodestotal; ++j)
1810  {
1811  InvR.SetValue(
1812  GetVertexMap(i), facemodearray[j],
1813  -R(GetVertexMap(i), facemodearray[j]));
1814  for (n = 0; n < nedgemodestotal; ++n)
1815  {
1816  MatrixValue = InvR.GetValue(GetVertexMap(i),
1817  facemodearray[j])
1818  + R(GetVertexMap(i), edgemodearray[n])
1819  * R(edgemodearray[n], facemodearray[j]);
1820  InvR.SetValue(GetVertexMap(i),
1821  facemodearray[j],
1822  MatrixValue);
1823  }
1824  }
1825  }
1826 
1827  // Edge-face contributions
1828  for (i = 0; i < nedgemodestotal; ++i)
1829  {
1830  for (j = 0; j < nfacemodestotal; ++j)
1831  {
1832  InvR.SetValue(
1833  edgemodearray[i], facemodearray[j],
1834  -R(edgemodearray[i], facemodearray[j]));
1835  }
1836  }
1837 
1838  for (i = 0; i < nCoeffs; ++i)
1839  {
1840  InvR.SetValue(i, i, 1.0);
1841  }
1842 
1843  return m_inversetransformationmatrix;
1844  }
1845 
1846  Array<OneD, unsigned int> Expansion3D::v_GetEdgeInverseBoundaryMap(
1847  int eid)
1848  {
1849  int n, j;
1850  int nEdgeCoeffs;
1851  int nBndCoeffs = NumBndryCoeffs();
1852 
1853  Array<OneD, unsigned int> bmap(nBndCoeffs);
1854  GetBoundaryMap(bmap);
1855 
1856  // Map from full system to statically condensed system (i.e reverse
1857  // GetBoundaryMap)
1858  map<int, int> invmap;
1859  for (j = 0; j < nBndCoeffs; ++j)
1860  {
1861  invmap[bmap[j]] = j;
1862  }
1863 
1864  // Number of interior edge coefficients
1865  nEdgeCoeffs = GetEdgeNcoeffs(eid) - 2;
1866 
1867  const SpatialDomains::Geometry3DSharedPtr &geom = GetGeom3D();
1868 
1869  Array<OneD, unsigned int> edgemaparray(nEdgeCoeffs);
1870  StdRegions::Orientation eOrient =
1871  geom->GetEorient(eid);
1872  Array<OneD, unsigned int> maparray =
1873  Array<OneD, unsigned int>(nEdgeCoeffs);
1874  Array<OneD, int> signarray =
1875  Array<OneD, int>(nEdgeCoeffs, 1);
1876 
1877  // maparray is the location of the edge within the matrix
1878  GetEdgeInteriorMap(eid, eOrient, maparray, signarray);
1879 
1880  for (n = 0; n < nEdgeCoeffs; ++n)
1881  {
1882  edgemaparray[n] = invmap[maparray[n]];
1883  }
1884 
1885  return edgemaparray;
1886  }
1887 
1888  Array<OneD, unsigned int>
1889  Expansion3D::v_GetFaceInverseBoundaryMap(
1890  int fid,
1891  StdRegions::Orientation faceOrient)
1892  {
1893  int n, j;
1894  int nFaceCoeffs;
1895 
1896  int nBndCoeffs = NumBndryCoeffs();
1897 
1898  Array<OneD, unsigned int> bmap(nBndCoeffs);
1899  GetBoundaryMap(bmap);
1900 
1901  // Map from full system to statically condensed system (i.e reverse
1902  // GetBoundaryMap)
1903  map<int, int> reversemap;
1904  for (j = 0; j < bmap.num_elements(); ++j)
1905  {
1906  reversemap[bmap[j]] = j;
1907  }
1908 
1909  // Number of interior face coefficients
1910  nFaceCoeffs = GetFaceIntNcoeffs(fid);
1911 
1912  Array<OneD, unsigned int> facemaparray(nFaceCoeffs);
1913  StdRegions::Orientation fOrient;
1914  Array<OneD, unsigned int> maparray =
1915  Array<OneD, unsigned int>(nFaceCoeffs);
1916  Array<OneD, int> signarray =
1917  Array<OneD, int>(nFaceCoeffs, 1);
1918 
1919  if(faceOrient == StdRegions::eNoOrientation)
1920  {
1921  fOrient = GetForient(fid);
1922  }
1923  else
1924  {
1925  fOrient = faceOrient;
1926  }
1927 
1928  // maparray is the location of the face within the matrix
1929  GetFaceInteriorMap(fid, fOrient, maparray, signarray);
1930 
1931  for (n = 0; n < nFaceCoeffs; ++n)
1932  {
1933  facemaparray[n] = reversemap[maparray[n]];
1934  }
1935 
1936  return facemaparray;
1937  }
1938 
1939  StdRegions::Orientation Expansion3D::v_GetForient(int face)
1940  {
1941  return m_geom->GetForient(face);
1942  }
1943 
1944  /**
1945  * \brief Returns the physical values at the quadrature points of a face
1946  * Wrapper function to v_GetFacePhysVals
1947  */
1948  void Expansion3D::v_GetTracePhysVals(
1949  const int face,
1950  const StdRegions::StdExpansionSharedPtr &FaceExp,
1951  const Array<OneD, const NekDouble> &inarray,
1952  Array<OneD, NekDouble> &outarray,
1953  StdRegions::Orientation orient)
1954  {
1955  v_GetFacePhysVals(face,FaceExp,inarray,outarray,orient);
1956  }
1957 
1958  void Expansion3D::v_GetFacePhysVals(
1959  const int face,
1960  const StdRegions::StdExpansionSharedPtr &FaceExp,
1961  const Array<OneD, const NekDouble> &inarray,
1962  Array<OneD, NekDouble> &outarray,
1963  StdRegions::Orientation orient)
1964  {
1965 
1966  if (orient == StdRegions::eNoOrientation)
1967  {
1968  orient = GetForient(face);
1969  }
1970 
1971  int nq0 = FaceExp->GetNumPoints(0);
1972  int nq1 = FaceExp->GetNumPoints(1);
1973 
1974  int nfacepts = GetFaceNumPoints(face);
1975  int dir0 = GetGeom3D()->GetDir(face,0);
1976  int dir1 = GetGeom3D()->GetDir(face,1);
1977 
1978  Array<OneD,NekDouble> o_tmp (nfacepts);
1979  Array<OneD,NekDouble> o_tmp2(FaceExp->GetTotPoints());
1980  Array<OneD, int> faceids;
1981 
1982  // Get local face pts and put into o_tmp
1983  GetFacePhysMap(face,faceids);
1984  Vmath::Gathr(faceids.num_elements(),inarray,faceids,o_tmp);
1985 
1986 
1987  int to_id0,to_id1;
1988 
1989  if(orient < StdRegions::eDir1FwdDir2_Dir2FwdDir1)
1990  {
1991  to_id0 = 0;
1992  to_id1 = 1;
1993  }
1994  else // transpose points key evaluation
1995  {
1996  to_id0 = 1;
1997  to_id1 = 0;
1998  }
1999 
2000  // interpolate to points distrbution given in FaceExp
2001  LibUtilities::Interp2D(m_base[dir0]->GetPointsKey(),
2002  m_base[dir1]->GetPointsKey(),
2003  o_tmp.get(),
2004  FaceExp->GetBasis(to_id0)->GetPointsKey(),
2005  FaceExp->GetBasis(to_id1)->GetPointsKey(),
2006  o_tmp2.get());
2007 
2008  // Reshuffule points as required and put into outarray.
2009  ReOrientFacePhysMap(FaceExp->GetNverts(),orient,nq0,nq1,faceids);
2010  Vmath::Scatr(nq0*nq1,o_tmp2,faceids,outarray);
2011  }
2012 
2013  void Expansion3D::ReOrientFacePhysMap(const int nvert,
2014  const StdRegions::Orientation orient,
2015  const int nq0, const int nq1,
2016  Array<OneD, int> &idmap)
2017  {
2018  if(nvert == 3)
2019  {
2020  ReOrientTriFacePhysMap(orient,nq0,nq1,idmap);
2021  }
2022  else
2023  {
2024  ReOrientQuadFacePhysMap(orient,nq0,nq1,idmap);
2025  }
2026  }
2027 
2028  void Expansion3D::ReOrientTriFacePhysMap(const StdRegions::Orientation orient,
2029  const int nq0,
2030  const int nq1,
2031  Array<OneD, int> &idmap)
2032  {
2033  if(idmap.num_elements() != nq0*nq1)
2034  {
2035  idmap = Array<OneD,int>(nq0*nq1);
2036  }
2037 
2038  switch(orient)
2039  {
2040  case StdRegions::eDir1FwdDir1_Dir2FwdDir2:
2041  // eseentially straight copy
2042  for(int i = 0; i < nq0*nq1; ++i)
2043  {
2044  idmap[i] = i;
2045  }
2046  break;
2047  case StdRegions::eDir1BwdDir1_Dir2FwdDir2:
2048  // reverse
2049  for (int j = 0; j < nq1; ++j)
2050  {
2051  for(int i = 0; i < nq0; ++i)
2052  {
2053  idmap[j*nq0+i] = nq0-1-i +j*nq0;
2054  }
2055  }
2056  break;
2057  default:
2058  ASSERTL0(false,"Case not supposed to be used in this function");
2059  }
2060  }
2061 
2062 
2063  void Expansion3D::ReOrientQuadFacePhysMap(const StdRegions::Orientation orient,
2064  const int nq0,
2065  const int nq1,
2066  Array<OneD, int> &idmap)
2067  {
2068  if(idmap.num_elements() != nq0*nq1)
2069  {
2070  idmap = Array<OneD,int>(nq0*nq1);
2071  }
2072 
2073 
2074  switch(orient)
2075  {
2076  case StdRegions::eDir1FwdDir1_Dir2FwdDir2:
2077  // eseentially straight copy
2078  for(int i = 0; i < nq0*nq1; ++i)
2079  {
2080  idmap[i] = i;
2081  }
2082  break;
2083  case StdRegions::eDir1BwdDir1_Dir2FwdDir2:
2084  {
2085  //Direction A negative and B positive
2086  for (int j = 0; j < nq1; j++)
2087  {
2088  for (int i =0; i < nq0; ++i)
2089  {
2090  idmap[j*nq0+i] = nq0-1-i + j*nq0;
2091  }
2092  }
2093  }
2094  break;
2095  case StdRegions::eDir1FwdDir1_Dir2BwdDir2:
2096  {
2097  //Direction A positive and B negative
2098  for (int j = 0; j < nq1; j++)
2099  {
2100  for (int i =0; i < nq0; ++i)
2101  {
2102  idmap[j*nq0+i] = nq0*(nq1-1-j)+i;
2103  }
2104  }
2105  }
2106  case StdRegions::eDir1BwdDir1_Dir2BwdDir2:
2107  {
2108  //Direction A negative and B negative
2109  for (int j = 0; j < nq1; j++)
2110  {
2111  for (int i =0; i < nq0; ++i)
2112  {
2113  idmap[j*nq0+i] = nq0*nq1-1-j*nq0 -i;
2114  }
2115  }
2116  }
2117  case StdRegions::eDir1FwdDir2_Dir2FwdDir1:
2118  {
2119  //Transposed, Direction A and B positive
2120  for (int i =0; i < nq0; ++i)
2121  {
2122  for (int j = 0; j < nq1; ++j)
2123  {
2124  idmap[i*nq1+j] = i + j*nq0;
2125  }
2126  }
2127  }
2128  break;
2129  case StdRegions::eDir1FwdDir2_Dir2BwdDir1:
2130  {
2131  //Transposed, Direction A positive and B negative
2132  for (int i =0; i < nq0; ++i)
2133  {
2134  for (int j = 0; j < nq1; ++j)
2135  {
2136  idmap[i*nq1+j] = nq0-1-i + j*nq0;
2137  }
2138  }
2139  }
2140  break;
2141  case StdRegions::eDir1BwdDir2_Dir2FwdDir1:
2142  {
2143  //Transposed, Direction A negative and B positive
2144  for (int i =0; i < nq0; ++i)
2145  {
2146  for (int j = 0; j < nq1; ++j)
2147  {
2148  idmap[i*nq1+j] = i+nq0*(nq1-1)-j*nq0;
2149  }
2150  }
2151  break;
2152  }
2153  case StdRegions::eDir1BwdDir2_Dir2BwdDir1:
2154  {
2155  //Transposed, Direction A and B negative
2156  for (int i =0; i < nq0; ++i)
2157  {
2158  for (int j = 0; j < nq1; ++j)
2159  {
2160  idmap[i*nq1+j] = nq0*nq1-1-i-j*nq0;
2161  }
2162  }
2163  }
2164  break;
2165  default:
2166  ASSERTL0(false,"Unknow orientation");
2167  break;
2168  }
2169  }
2170 
2171  void Expansion3D::v_NormVectorIProductWRTBase(
2172  const Array<OneD, const Array<OneD, NekDouble> > &Fvec,
2173  Array<OneD, NekDouble> &outarray)
2174  {
2175  NormVectorIProductWRTBase(Fvec[0], Fvec[1], Fvec[2], outarray);
2176  }
2177  } //end of namespace
2178 } //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:198
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:191
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:644
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:27
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:485
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:442
Nektar::LibUtilities::eModified_A
Principle Modified Functions .
Definition: BasisType.h:49
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:227
Vmath::Scatr
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition: Vmath.cpp:673
Nektar::Transpose
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Definition: BlockMatrix.hpp:219
Expansion3D.h
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:396
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:250
InterpCoeff.h
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
Nektar::StdRegions::StdExpansionSharedPtr
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1878
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
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:183
Nektar::LocalRegions::Expansion2DSharedPtr
boost::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:49
Nektar::StdRegions::NullConstFactorMap
static ConstFactorMap NullConstFactorMap
Definition: StdRegions.hpp:253
Nektar::SpatialDomains::Geometry3DSharedPtr
boost::shared_ptr< Geometry3D > Geometry3DSharedPtr
Definition: Geometry3D.h:52
Generated on Wed Mar 8 2017 12:28:54 for Nektar++ by   doxygen 1.8.8