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