Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
StdExpansion.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File StdExpansion.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: Definition of methods in class StdExpansion which is
33 // the base class to all expansion shapes
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
37 #include <StdRegions/StdExpansion.h>
38 #include <LibUtilities/Foundations/ManagerAccess.h> // for BasisManager, etc
39 
40 namespace Nektar
41 {
42  namespace StdRegions
43  {
44  StdExpansion::StdExpansion(void):
45  m_elmt_id(0),
46  m_ncoeffs(0)
47  {
48  }
49 
50  StdExpansion::StdExpansion(const int numcoeffs, const int numbases,
51  const LibUtilities::BasisKey &Ba,
52  const LibUtilities::BasisKey &Bb,
53  const LibUtilities::BasisKey &Bc):
54  m_base(numbases),
55  m_elmt_id(0),
56  m_ncoeffs(numcoeffs),
57  m_stdMatrixManager(
58  boost::bind(&StdExpansion::CreateStdMatrix, this, _1),
59  std::string("StdExpansionStdMatrix")),
60  m_stdStaticCondMatrixManager(
61  boost::bind(&StdExpansion::CreateStdStaticCondMatrix, this, _1),
62  std::string("StdExpansionStdStaticCondMatrix")),
63  m_IndexMapManager(
64  boost::bind(&StdExpansion::CreateIndexMap,this, _1),
65  std::string("StdExpansionIndexMap"))
66  {
67  switch(m_base.num_elements())
68  {
69  case 3:
70  ASSERTL2(Bc!=LibUtilities::NullBasisKey,
71  "NULL Basis attempting to be used.");
72  m_base[2] = LibUtilities::BasisManager()[Bc];
73 
74  case 2:
75  ASSERTL2(Bb!=LibUtilities::NullBasisKey,
76  "NULL Basis attempting to be used.");
77 
78  m_base[1] = LibUtilities::BasisManager()[Bb];
79  case 1:
80  ASSERTL2(Ba!=LibUtilities::NullBasisKey,
81  "NULL Basis attempting to be used.");
82  m_base[0] = LibUtilities::BasisManager()[Ba];
83  break;
84  default:
85  break;
86 // ASSERTL0(false, "numbases incorrectly specified");
87  };
88 
89  } //end constructor
90 
91 
92  StdExpansion::StdExpansion(const StdExpansion &T):
93  m_base(T.m_base),
94  m_elmt_id(T.m_elmt_id),
95  m_ncoeffs(T.m_ncoeffs),
96  m_stdMatrixManager(T.m_stdMatrixManager),
97  m_stdStaticCondMatrixManager(T.m_stdStaticCondMatrixManager),
98  m_IndexMapManager(T.m_IndexMapManager)
99  {
100  }
101 
102  StdExpansion::~StdExpansion()
103  {
104  }
105 
106  NekDouble StdExpansion::Linf(const Array<OneD, const NekDouble>& phys,
107  const Array<OneD, const NekDouble>& sol)
108  {
109  NekDouble val;
110  int ntot = GetTotPoints();
111  Array<OneD, NekDouble> wsp(ntot);
112 
113  if(sol == NullNekDouble1DArray)
114  {
115  Vmath::Vabs(ntot, phys, 1, wsp, 1);
116  }
117  else
118  {
119  Vmath::Vsub(ntot, sol, 1, phys, 1, wsp, 1);
120  Vmath::Vabs(ntot, wsp, 1, wsp, 1);
121  }
122 
123  val = Vmath::Vamax(ntot, wsp, 1);
124 
125  return val;
126  }
127 
128  NekDouble StdExpansion::L2(const Array<OneD, const NekDouble>& phys,
129  const Array<OneD, const NekDouble>& sol)
130  {
131  NekDouble val;
132  int ntot = GetTotPoints();
133  Array<OneD, NekDouble> wsp(ntot);
134 
135  if (sol.num_elements() == 0)
136  {
137  Vmath::Vmul(ntot, phys, 1, phys, 1, wsp, 1);
138  }
139  else
140  {
141  Vmath::Vsub(ntot, sol, 1, phys, 1, wsp, 1);
142  Vmath::Vmul(ntot, wsp, 1, wsp, 1, wsp, 1);
143  }
144 
145  val = v_Integral(wsp);
146 
147  // if val too small, sqrt returns nan.
148  if (fabs(val) < NekConstants::kNekSqrtTol*NekConstants::kNekSqrtTol)
149  {
150  return 0.0;
151  }
152  else
153  {
154  return sqrt(val);
155  }
156  }
157 
158  NekDouble StdExpansion::H1(const Array<OneD, const NekDouble>& phys,
159  const Array<OneD, const NekDouble>& sol)
160  {
161  int i;
162  NekDouble val;
163  int ntot = GetTotPoints();
164  int coordim = v_GetCoordim();
165  Array<OneD, NekDouble> wsp(3*ntot);
166  Array<OneD, NekDouble> wsp_deriv = wsp + ntot;
167  Array<OneD, NekDouble> sum = wsp_deriv + ntot;
168 
169  if(sol == NullNekDouble1DArray)
170  {
171  Vmath::Vcopy(ntot,phys, 1, wsp, 1);
172  Vmath::Vmul(ntot, phys, 1, phys, 1, sum, 1);
173  }
174  else
175  {
176  Vmath::Vsub(ntot, sol, 1, phys, 1, wsp, 1);
177  Vmath::Vmul(ntot, wsp, 1, wsp, 1, sum, 1);
178  }
179 
180 
181  for(i = 0; i < coordim; ++i)
182  {
183  v_PhysDeriv(i,wsp,wsp_deriv);
184  Vmath::Vvtvp(ntot,wsp_deriv,1,wsp_deriv,1,sum,1,sum,1);
185  }
186 
187  val = sqrt(v_Integral(sum));
188 
189  return val;
190  }
191 
192 
193  DNekBlkMatSharedPtr StdExpansion::CreateStdStaticCondMatrix(const StdMatrixKey &mkey)
194  {
195  DNekBlkMatSharedPtr returnval;
196 
197  DNekMatSharedPtr mat = GetStdMatrix(mkey);
198  int nbdry = NumBndryCoeffs(); // also checks to see if this is a boundary interior decomposed expansion
199  int nint = m_ncoeffs - nbdry;
200  DNekMatSharedPtr A = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nbdry);
201  DNekMatSharedPtr B = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nint);
202  DNekMatSharedPtr C = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nbdry);
203  DNekMatSharedPtr D = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nint);
204 
205  int i,j;
206 
207  Array<OneD,unsigned int> bmap(nbdry);
208  Array<OneD,unsigned int> imap(nint);
209  GetBoundaryMap(bmap);
210  GetInteriorMap(imap);
211 
212  for(i = 0; i < nbdry; ++i)
213  {
214  for(j = 0; j < nbdry; ++j)
215  {
216  (*A)(i,j) = (*mat)(bmap[i],bmap[j]);
217  }
218 
219  for(j = 0; j < nint; ++j)
220  {
221  (*B)(i,j) = (*mat)(bmap[i],imap[j]);
222  }
223  }
224 
225  for(i = 0; i < nint; ++i)
226  {
227  for(j = 0; j < nbdry; ++j)
228  {
229  (*C)(i,j) = (*mat)(imap[i],bmap[j]);
230  }
231 
232  for(j = 0; j < nint; ++j)
233  {
234  (*D)(i,j) = (*mat)(imap[i],imap[j]);
235  }
236  }
237 
238  // Calculate static condensed system
239  if(nint)
240  {
241  D->Invert();
242  (*B) = (*B)*(*D);
243  (*A) = (*A) - (*B)*(*C);
244  }
245 
246  // set up block matrix system
247  Array<OneD, unsigned int> exp_size(2);
248  exp_size[0] = nbdry;
249  exp_size[1] = nint;
250  returnval = MemoryManager<DNekBlkMat>::AllocateSharedPtr(exp_size,exp_size);
251 
252  returnval->SetBlock(0,0,A);
253  returnval->SetBlock(0,1,B);
254  returnval->SetBlock(1,0,C);
255  returnval->SetBlock(1,1,D);
256 
257  return returnval;
258  }
259 
260  IndexMapValuesSharedPtr StdExpansion::CreateIndexMap(const IndexMapKey &ikey)
261  {
262  IndexMapValuesSharedPtr returnval;
263 
264  IndexMapType itype = ikey.GetIndexMapType();
265 
266  int entity = ikey.GetIndexEntity();
267 
268  Orientation orient = ikey.GetIndexOrientation();
269 
270  Array<OneD,unsigned int> map;
271  Array<OneD,int> sign;
272 
273  switch(itype)
274  {
275  case eEdgeToElement:
276  {
277  v_GetEdgeToElementMap(entity,orient,map,sign);
278  }
279  break;
280  case eFaceToElement:
281  {
282  v_GetFaceToElementMap(entity,orient,map,sign);
283  }
284  break;
285  case eEdgeInterior:
286  {
287  v_GetEdgeInteriorMap(entity,orient,map,sign);
288  }
289  break;
290  case eFaceInterior:
291  {
292  v_GetFaceInteriorMap(entity,orient,map,sign);
293  }
294  break;
295  case eBoundary:
296  {
297  ASSERTL0(false,"Boundary Index Map not implemented yet.");
298  }
299  break;
300  case eVertex:
301  {
302  ASSERTL0(false,"Vertex Index Map not implemented yet.");
303  }
304  break;
305  default:
306  {
307  ASSERTL0(false,"The Index Map you are requiring is not between the possible options.");
308  }
309  }
310 
311  returnval = MemoryManager<IndexMapValues>::AllocateSharedPtr(map.num_elements());
312 
313  for(int i = 0; i < map.num_elements(); i++)
314  {
315  (*returnval)[i].index = map[i];
316  (*returnval)[i].sign = sign[i];
317  }
318 
319  return returnval;
320  }
321 
322  DNekMatSharedPtr StdExpansion::CreateGeneralMatrix(const StdMatrixKey &mkey)
323  {
324  int i;
325  DNekMatSharedPtr returnval;
326 
327  switch(mkey.GetMatrixType())
328  {
329  case eInvMass:
330  {
331  StdMatrixKey masskey(eMass,mkey.GetShapeType(),*this,NullConstFactorMap,NullVarCoeffMap,mkey.GetNodalPointsType());
332  DNekMatSharedPtr mmat = GetStdMatrix(masskey);
333 
334  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(*mmat); //Populate standard mass matrix.
335  returnval->Invert();
336  }
337  break;
338  case eInvNBasisTrans:
339  {
340  StdMatrixKey tmpkey(eNBasisTrans,mkey.GetShapeType(),*this,NullConstFactorMap,NullVarCoeffMap,mkey.GetNodalPointsType());
341  DNekMatSharedPtr tmpmat = GetStdMatrix(tmpkey);
342  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(*tmpmat); //Populate matrix.
343  returnval->Invert();
344  }
345  break;
346  case eBwdTrans:
347  {
348  int nq = GetTotPoints();
349  Array<OneD, NekDouble> tmpin(m_ncoeffs);
350  Array<OneD, NekDouble> tmpout(nq);
351 
352  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(nq,m_ncoeffs);
353 
354  for(int i=0; i<m_ncoeffs; ++i)
355  {
356  Vmath::Zero(m_ncoeffs, tmpin, 1);
357  tmpin[i] = 1.0;
358 
359  BwdTrans_SumFac(tmpin,tmpout);
360 
361  Vmath::Vcopy(nq,tmpout.get(),1,
362  returnval->GetRawPtr()+i*nq,1);
363  }
364  }
365  break;
366  case eIProductWRTBase:
367  {
368  int nq = GetTotPoints();
369  Array<OneD, NekDouble> tmpin(nq);
370  Array<OneD, NekDouble> tmpout(m_ncoeffs);
371 
372  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(m_ncoeffs,nq);
373 
374  for(i=0; i < nq; ++i)
375  {
376  Vmath::Zero(nq, tmpin, 1);
377  tmpin[i] = 1.0;
378 
379  IProductWRTBase_SumFac(tmpin,tmpout);
380 
381  Vmath::Vcopy(m_ncoeffs,tmpout.get(),1,
382  returnval->GetRawPtr()+i*m_ncoeffs,1);
383  }
384  }
385  break;
386  case eIProductWRTDerivBase0:
387  {
388  int nq = GetTotPoints();
389  Array<OneD, NekDouble> tmpin(nq);
390  Array<OneD, NekDouble> tmpout(m_ncoeffs);
391 
392  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(nq,m_ncoeffs);
393 
394  for(i=0; i < nq; ++i)
395  {
396  Vmath::Zero(nq, tmpin, 1);
397  tmpin[i] = 1.0;
398 
399  IProductWRTDerivBase_SumFac(0,tmpin,tmpout);
400 
401  Vmath::Vcopy(m_ncoeffs,tmpout.get(),1,
402  returnval->GetRawPtr()+i*m_ncoeffs,1);
403  }
404  }
405  break;
406  case eIProductWRTDerivBase1:
407  {
408  int nq = GetTotPoints();
409  Array<OneD, NekDouble> tmpin(nq);
410  Array<OneD, NekDouble> tmpout(m_ncoeffs);
411 
412  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(nq,m_ncoeffs);
413 
414  for(i=0; i < nq; ++i)
415  {
416  Vmath::Zero(nq, tmpin, 1);
417  tmpin[i] = 1.0;
418 
419  IProductWRTDerivBase_SumFac(1,tmpin,tmpout);
420 
421  Vmath::Vcopy(m_ncoeffs,tmpout.get(),1,
422  returnval->GetRawPtr()+i*m_ncoeffs,1);
423  }
424  }
425  break;
426  case eIProductWRTDerivBase2:
427  {
428  int nq = GetTotPoints();
429  Array<OneD, NekDouble> tmpin(nq);
430  Array<OneD, NekDouble> tmpout(m_ncoeffs);
431 
432  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(nq,m_ncoeffs);
433 
434  for(i=0; i < nq; ++i)
435  {
436  Vmath::Zero(nq, tmpin, 1);
437  tmpin[i] = 1.0;
438 
439  IProductWRTDerivBase_SumFac(2,tmpin,tmpout);
440 
441  Vmath::Vcopy(m_ncoeffs,tmpout.get(),1,
442  returnval->GetRawPtr()+i*m_ncoeffs,1);
443  }
444  }
445  break;
446  case eMass:
447  case eHelmholtz:
448  case eLaplacian:
449  case eLaplacian00:
450  case eLaplacian01:
451  case eLaplacian02:
452  case eLaplacian11:
453  case eLaplacian12:
454  case eLaplacian22:
455  case eWeakDeriv0:
456  case eWeakDeriv1:
457  case eWeakDeriv2:
458  case eWeakDirectionalDeriv:
459  case eMassLevelCurvature:
460  case eLinearAdvectionReaction:
461  case eLinearAdvectionDiffusionReaction:
462  {
463  Array<OneD, NekDouble> tmp(m_ncoeffs);
464  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(m_ncoeffs,m_ncoeffs);
465  DNekMat &Mat = *returnval;
466 
467  for(i=0; i < m_ncoeffs; ++i)
468  {
469  Vmath::Zero(m_ncoeffs, tmp, 1);
470  tmp[i] = 1.0;
471 
472  GeneralMatrixOp_MatFree(tmp,tmp,mkey);
473 
474  Vmath::Vcopy(m_ncoeffs,&tmp[0],1,
475  &(Mat.GetPtr())[0]+i*m_ncoeffs,1);
476  }
477  }
478  break;
479  default:
480  {
481  NEKERROR(ErrorUtil::efatal, "This type of matrix can not be created using a general approach");
482  }
483  break;
484  }
485 
486  return returnval;
487  }
488 
489  void StdExpansion::GeneralMatrixOp(const Array<OneD, const NekDouble> &inarray,
490  Array<OneD,NekDouble> &outarray,
491  const StdMatrixKey &mkey)
492  {
493  switch(mkey.GetMatrixType())
494  {
495  case eMass:
496  MassMatrixOp(inarray,outarray,mkey);
497  break;
498  case eWeakDeriv0:
499  WeakDerivMatrixOp(0,inarray,outarray,mkey);
500  break;
501  case eWeakDeriv1:
502  WeakDerivMatrixOp(1,inarray,outarray,mkey);
503  break;
504  case eWeakDeriv2:
505  WeakDerivMatrixOp(2,inarray,outarray,mkey);
506  break;
507  case eWeakDirectionalDeriv:
508  WeakDirectionalDerivMatrixOp(inarray,outarray,mkey);
509  break;
510  case eMassLevelCurvature:
511  MassLevelCurvatureMatrixOp(inarray,outarray,mkey);
512  break;
513  case eLinearAdvectionReaction:
514  LinearAdvectionDiffusionReactionMatrixOp(inarray,outarray,mkey,false);
515  break;
516  case eLinearAdvectionDiffusionReaction:
517  LinearAdvectionDiffusionReactionMatrixOp(inarray,outarray,mkey);
518  break;
519  case eLaplacian:
520  LaplacianMatrixOp(inarray,outarray,mkey);
521  break;
522  case eLaplacian00:
523  LaplacianMatrixOp(0,0,inarray,outarray,mkey);
524  break;
525  case eLaplacian01:
526  LaplacianMatrixOp(0,1,inarray,outarray,mkey);
527  break;
528  case eLaplacian02:
529  LaplacianMatrixOp(0,2,inarray,outarray,mkey);
530  break;
531  case eLaplacian10:
532  LaplacianMatrixOp(1,0,inarray,outarray,mkey);
533  break;
534  case eLaplacian11:
535  LaplacianMatrixOp(1,1,inarray,outarray,mkey);
536  break;
537  case eLaplacian12:
538  LaplacianMatrixOp(1,2,inarray,outarray,mkey);
539  break;
540  case eLaplacian20:
541  LaplacianMatrixOp(2,0,inarray,outarray,mkey);
542  break;
543  case eLaplacian21:
544  LaplacianMatrixOp(2,1,inarray,outarray,mkey);
545  break;
546  case eLaplacian22:
547  LaplacianMatrixOp(2,2,inarray,outarray,mkey);
548  break;
549  case eHelmholtz:
550  HelmholtzMatrixOp(inarray,outarray,mkey);
551  break;
552  default:
553  NEKERROR(ErrorUtil::efatal, "This matrix does not have an operator");
554  break;
555  }
556  }
557 
558  void StdExpansion::GeneralMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
559  Array<OneD,NekDouble> &outarray,
560  const StdMatrixKey &mkey)
561  {
562  switch(mkey.GetMatrixType())
563  {
564  case eMass:
565  MassMatrixOp_MatFree(inarray,outarray,mkey);
566  break;
567  case eWeakDeriv0:
568  WeakDerivMatrixOp_MatFree(0,inarray,outarray,mkey);
569  break;
570  case eWeakDeriv1:
571  WeakDerivMatrixOp_MatFree(1,inarray,outarray,mkey);
572  break;
573  case eWeakDeriv2:
574  WeakDerivMatrixOp_MatFree(2,inarray,outarray,mkey);
575  break;
576  case eWeakDirectionalDeriv:
577  WeakDirectionalDerivMatrixOp_MatFree(inarray,outarray,mkey);
578  break;
579  case eMassLevelCurvature:
580  MassLevelCurvatureMatrixOp_MatFree(inarray,outarray,mkey);
581  break;
582  case eLinearAdvectionReaction:
583  LinearAdvectionDiffusionReactionMatrixOp_MatFree(inarray,outarray,mkey,false);
584  break;
585  case eLinearAdvectionDiffusionReaction:
586  LinearAdvectionDiffusionReactionMatrixOp_MatFree(inarray,outarray,mkey);
587  break;
588  case eLaplacian:
589  LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
590  break;
591  case eLaplacian00:
592  LaplacianMatrixOp_MatFree(0,0,inarray,outarray,mkey);
593  break;
594  case eLaplacian01:
595  LaplacianMatrixOp_MatFree(0,1,inarray,outarray,mkey);
596  break;
597  case eLaplacian02:
598  LaplacianMatrixOp_MatFree(0,2,inarray,outarray,mkey);
599  break;
600  case eLaplacian10:
601  LaplacianMatrixOp_MatFree(1,0,inarray,outarray,mkey);
602  break;
603  case eLaplacian11:
604  LaplacianMatrixOp_MatFree(1,1,inarray,outarray,mkey);
605  break;
606  case eLaplacian12:
607  LaplacianMatrixOp_MatFree(1,2,inarray,outarray,mkey);
608  break;
609  case eLaplacian20:
610  LaplacianMatrixOp_MatFree(2,0,inarray,outarray,mkey);
611  break;
612  case eLaplacian21:
613  LaplacianMatrixOp_MatFree(2,1,inarray,outarray,mkey);
614  break;
615  case eLaplacian22:
616  LaplacianMatrixOp_MatFree(2,2,inarray,outarray,mkey);
617  break;
618  case eHelmholtz:
619  HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
620  break;
621  default:
622  NEKERROR(ErrorUtil::efatal, "This matrix does not have an operator");
623  break;
624  }
625  }
626 
627  void StdExpansion::MassMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
628  Array<OneD,NekDouble> &outarray,
629  const StdMatrixKey &mkey)
630  {
631  int nq = GetTotPoints();
632  Array<OneD, NekDouble> tmp(nq);
633 
634  v_BwdTrans(inarray,tmp);
635 
636  if(mkey.HasVarCoeff(eVarCoeffMass))
637  {
638  Vmath::Vmul(nq, mkey.GetVarCoeff(eVarCoeffMass), 1, tmp, 1, tmp, 1);
639  }
640 
641  v_IProductWRTBase(tmp, outarray);
642  }
643 
644  void StdExpansion::LaplacianMatrixOp_MatFree(const int k1, const int k2,
645  const Array<OneD, const NekDouble> &inarray,
646  Array<OneD,NekDouble> &outarray,
647  const StdMatrixKey &mkey)
648  {
649  ASSERTL1(k1 >= 0 && k1 < GetCoordim(),"invalid first argument");
650  ASSERTL1(k2 >= 0 && k2 < GetCoordim(),"invalid second argument");
651 
652  int nq = GetTotPoints();
653  Array<OneD, NekDouble> tmp(nq);
654  Array<OneD, NekDouble> dtmp(nq);
655  VarCoeffType varcoefftypes[3][3]
656  = { {eVarCoeffD00, eVarCoeffD01, eVarCoeffD02},
657  {eVarCoeffD01, eVarCoeffD11, eVarCoeffD12},
658  {eVarCoeffD02, eVarCoeffD12, eVarCoeffD22}
659  };
660 
661  v_BwdTrans(inarray,tmp);
662  v_PhysDeriv(k2,tmp,dtmp);
663  if (mkey.GetNVarCoeff())
664  {
665  if (k1 == k2)
666  {
667  // By default, k1 == k2 has \sigma = 1 (diagonal entries)
668  if(mkey.HasVarCoeff(varcoefftypes[k1][k1]))
669  {
670  Vmath::Vmul(nq, mkey.GetVarCoeff(varcoefftypes[k1][k1]), 1, dtmp, 1, dtmp, 1);
671  }
672  v_IProductWRTDerivBase(k1, dtmp, outarray);
673  }
674  else
675  {
676  // By default, k1 != k2 has \sigma = 0 (off-diagonal entries)
677  if(mkey.HasVarCoeff(varcoefftypes[k1][k2]))
678  {
679  Vmath::Vmul(nq, mkey.GetVarCoeff(varcoefftypes[k1][k2]), 1, dtmp, 1, dtmp, 1);
680  v_IProductWRTDerivBase(k1, dtmp, outarray);
681  }
682  else
683  {
684  Vmath::Zero(GetNcoeffs(), outarray, 1);
685  }
686  }
687  }
688  else
689  {
690  // Multiply by svv tensor
691  if(mkey.ConstFactorExists(eFactorSVVCutoffRatio))
692  {
693  SVVLaplacianFilter(dtmp,mkey);
694  }
695  v_IProductWRTDerivBase(k1, dtmp, outarray);
696  }
697  }
698 
699  void StdExpansion::LaplacianMatrixOp_MatFree_GenericImpl(const Array<OneD, const NekDouble> &inarray,
700  Array<OneD,NekDouble> &outarray,
701  const StdMatrixKey &mkey)
702  {
703  const int dim = GetCoordim();
704 
705  int i,j;
706 
707  Array<OneD,NekDouble> store(m_ncoeffs);
708  Array<OneD,NekDouble> store2(m_ncoeffs,0.0);
709 
710  if(mkey.GetNVarCoeff() == 0)
711  {
712  // just call diagonal matrix form of laplcian operator
713  for(i = 0; i < dim; ++i)
714  {
715  LaplacianMatrixOp(i,i,inarray,store,mkey);
716  Vmath::Vadd(m_ncoeffs, store, 1, store2, 1, store2, 1);
717  }
718  }
719  else
720  {
721  const MatrixType mtype[3][3]
722  = {{eLaplacian00,eLaplacian01,eLaplacian02},
723  {eLaplacian01,eLaplacian11,eLaplacian12},
724  {eLaplacian02,eLaplacian12,eLaplacian22}};
725  StdMatrixKeySharedPtr mkeyij;
726 
727  for(i = 0; i < dim; i++)
728  {
729  for(j = 0; j < dim; j++)
730  {
731  mkeyij = MemoryManager<StdMatrixKey>::AllocateSharedPtr(mkey,mtype[i][j]);
732  LaplacianMatrixOp(i,j,inarray,store,*mkeyij);
733  Vmath::Vadd(m_ncoeffs, store, 1, store2, 1, store2, 1);
734  }
735  }
736  }
737 
738  Vmath::Vcopy(m_ncoeffs,store2.get(),1,outarray.get(),1);
739  }
740 
741  void StdExpansion::WeakDerivMatrixOp_MatFree(const int k1,
742  const Array<OneD, const NekDouble> &inarray,
743  Array<OneD,NekDouble> &outarray,
744  const StdMatrixKey &mkey)
745  {
746  Array<OneD, NekDouble> tmp(GetTotPoints());
747  int nq = GetTotPoints();
748 
749  v_BwdTrans(inarray,tmp);
750  v_PhysDeriv(k1,tmp,tmp);
751 
752  VarCoeffType keys[] = {eVarCoeffD00, eVarCoeffD11, eVarCoeffD22};
753  if(mkey.HasVarCoeff(keys[k1]))
754  {
755  Vmath::Vmul(nq, &(mkey.GetVarCoeff(keys[k1]))[0], 1, &tmp[0], 1, &tmp[0], 1);
756  }
757 
758  v_IProductWRTBase(tmp, outarray);
759  }
760 
761  void StdExpansion::WeakDirectionalDerivMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
762  Array<OneD,NekDouble> &outarray,
763  const StdMatrixKey &mkey)
764  {
765  int nq = GetTotPoints();
766  // int varsize = ((mkey.GetVariableCoefficient(0)).num_elements())/dim;
767  Array<OneD, NekDouble> tmp(nq);
768 
769  v_BwdTrans(inarray,tmp);
770  // For Deformed mesh ==============
771  // if (varsize==nq)
772  // {
773  // v_PhysDirectionalDeriv(tmp,mkey.GetVariableCoefficient(0),tmp);
774  // }
775  //
776  // // For Regular mesh ==========
777  // else
778  // {
779  // ASSERTL0(false, "Wrong route");
780  // }
781 
782  v_IProductWRTBase(tmp, outarray);
783  }
784 
785  void StdExpansion::MassLevelCurvatureMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
786  Array<OneD,NekDouble> &outarray,
787  const StdMatrixKey &mkey)
788  {
789  ///@todo fix this
790  // int nqtot = GetTotPoints();
791  // int matrixid = mkey.GetMatrixID();
792  //
793  // NekDouble checkweight=0.0;
794  // Array<OneD, NekDouble> tmp(nqtot), tan(nqtot), dtan0(nqtot), dtan1(nqtot), weight(nqtot,0.0);
795  //
796  // int gmatnumber = (mkey.GetVariableCoefficient(1)).num_elements();
797  //
798  // v_BwdTrans(inarray,tmp);
799  //
800  // // weight = \grad \cdot tanvec
801  // for(int k = 0; k < GetCoordim(); ++k)
802  // {
803  // Vmath::Vcopy(nqtot, &(mkey.GetVariableCoefficient(0))[k*nqtot],
804  // 1, &tan[0], 1);
805  //
806  // // For Regular mesh ...
807  // if(gmatnumber==1)
808  // {
809  // // D_{/xi} and D_{/eta}
810  // v_PhysDeriv(0,tan,dtan0);
811  // v_PhysDeriv(1,tan,dtan1);
812  //
813  // // d v / d x_i = (d \xi / d x_i)*( d v / d \xi ) + (d \eta / d x_i)*( d v / d \eta )
814  // Vmath::Svtvp(nqtot,(mkey.GetVariableCoefficient(2*k+1))[0],&dtan0[0],1,&weight[0],1,&weight[0],1);
815  // Vmath::Svtvp(nqtot,(mkey.GetVariableCoefficient(2*k+2))[0],&dtan1[0],1,&weight[0],1,&weight[0],1);
816  // }
817  //
818  // // For Curved mesh ...
819  // else if(gmatnumber==nqtot)
820  // {
821  // // D_{x} and D_{y}
822  // v_PhysDeriv(k,tan,dtan0);
823  // Vmath::Vadd(nqtot,&dtan0[0],1,&weight[0],1,&weight[0],1);
824  // }
825  //
826  // else
827  // {
828  // ASSERTL1( ((gmatnumber=1) || (gmatnumber==nqtot) ), "Gmat is not in a right size");
829  // }
830  // }
831  //
832  // Vmath::Vmul(nqtot, &weight[0], 1, &tmp[0], 1, &tmp[0], 1);
833  // v_IProductWRTBase(tmp, outarray);
834  }
835 
836  void StdExpansion::LinearAdvectionDiffusionReactionMatrixOp_MatFree( const Array<OneD, const NekDouble> &inarray,
837  Array<OneD,NekDouble> &outarray,
838  const StdMatrixKey &mkey,
839  bool addDiffusionTerm)
840  {
841 
842  int i;
843  int ndir = mkey.GetNVarCoeff(); // assume num.r consts corresponds to directions
844  ASSERTL0(ndir,"Must define at least one advection velocity");
845 
846  NekDouble lambda = mkey.GetConstFactor(eFactorLambda);
847  int totpts = GetTotPoints();
848  Array<OneD, NekDouble> tmp(3*totpts);
849  Array<OneD, NekDouble> tmp_deriv = tmp + totpts;
850  Array<OneD, NekDouble> tmp_adv = tmp_deriv + totpts;
851 
852 
853  ASSERTL1(ndir <= GetCoordim(),"Number of constants is larger than coordinate dimensions");
854 
855  v_BwdTrans(inarray,tmp);
856 
857  VarCoeffType varcoefftypes[] = {eVarCoeffVelX, eVarCoeffVelY};
858 
859  //calculate u dx + v dy + ..
860  Vmath::Zero(totpts,tmp_adv,1);
861  for(i = 0; i < ndir; ++i)
862  {
863  v_PhysDeriv(i,tmp,tmp_deriv);
864  Vmath::Vvtvp(totpts,mkey.GetVarCoeff(varcoefftypes[i]),1,tmp_deriv,1,tmp_adv,1,tmp_adv,1);
865  }
866 
867  if(lambda) // add -lambda*u
868  {
869  Vmath::Svtvp(totpts,-lambda,tmp,1,tmp_adv,1,tmp_adv,1);
870  }
871 
872 
873  if(addDiffusionTerm)
874  {
875  Array<OneD, NekDouble> lap(m_ncoeffs);
876  StdMatrixKey mkeylap(eLaplacian,DetShapeType(),*this);
877  LaplacianMatrixOp(inarray,lap,mkeylap);
878 
879  v_IProductWRTBase(tmp_adv, outarray);
880  // Lap v - u.grad v + lambda*u
881  // => (grad u, grad v) + u.grad v - lambda*u
882  Vmath::Vadd(m_ncoeffs,lap,1,outarray,1,outarray,1);
883  }
884  else
885  {
886  v_IProductWRTBase(tmp_adv, outarray);
887  }
888 
889  }
890 
891 
892  void StdExpansion::HelmholtzMatrixOp_MatFree_GenericImpl(const Array<OneD, const NekDouble> &inarray,
893  Array<OneD,NekDouble> &outarray,
894  const StdMatrixKey &mkey)
895  {
896  NekDouble lambda = mkey.GetConstFactor(eFactorLambda);
897  Array<OneD,NekDouble> tmp(m_ncoeffs);
898  StdMatrixKey mkeymass(eMass,DetShapeType(),*this);
899  StdMatrixKey mkeylap(eLaplacian,DetShapeType(),*this);
900 
901  MassMatrixOp(inarray,tmp,mkeymass);
902  LaplacianMatrixOp(inarray,outarray,mkeylap);
903 
904  Blas::Daxpy(m_ncoeffs, lambda, tmp, 1, outarray, 1);
905  }
906 
907  void StdExpansion::BwdTrans_MatOp(const Array<OneD, const NekDouble>& inarray,
908  Array<OneD, NekDouble> &outarray)
909  {
910  int nq = GetTotPoints();
911  StdMatrixKey bwdtransmatkey(eBwdTrans,DetShapeType(),*this);
912  DNekMatSharedPtr bwdtransmat = GetStdMatrix(bwdtransmatkey);
913 
914  Blas::Dgemv('N',nq,m_ncoeffs,1.0,bwdtransmat->GetPtr().get(),
915  nq, inarray.get(), 1, 0.0, outarray.get(), 1);
916  }
917 
918  // VIRTUAL INLINE FUNCTIONS FROM HEADER FILE
919  void StdExpansion::SetUpPhysNormals(const int edge)
920  {
921  v_SetUpPhysNormals(edge);
922  }
923 
924  NekDouble StdExpansion::StdPhysEvaluate(const Array<OneD, const NekDouble> &Lcoord,
925  const Array<OneD, const NekDouble> &physvals)
926  {
927  return v_StdPhysEvaluate(Lcoord,physvals);
928  }
929 
930  int StdExpansion::v_GetElmtId(void)
931  {
932  return m_elmt_id;
933  }
934 
935  const Array<OneD, const NekDouble>& StdExpansion::v_GetPhysNormals(void)
936  {
937  NEKERROR(ErrorUtil::efatal, "This function is not valid for this class");
938  return NullNekDouble1DArray;
939  }
940 
941 
942  void StdExpansion::v_SetPhysNormals(Array<OneD, const NekDouble> &normal)
943  {
944  NEKERROR(ErrorUtil::efatal, "This function is not valid for this class");
945  }
946 
947  void StdExpansion::v_SetUpPhysNormals(const int edge)
948  {
949  NEKERROR(ErrorUtil::efatal, "This function is not valid for this class");
950  }
951 
952  int StdExpansion::v_CalcNumberOfCoefficients(const std::vector<unsigned int> &nummodes, int &modes_offset)
953  {
954  NEKERROR(ErrorUtil::efatal, "This function is not defined for this class");
955  return 0;
956  }
957 
958  void StdExpansion::v_ExtractDataToCoeffs(const NekDouble *data,
959  const std::vector<unsigned int > &nummodes,
960  const int nmode_offset,
961  NekDouble *coeffs)
962  {
963  NEKERROR(ErrorUtil::efatal, "This function is not defined for this class");
964  }
965 
966 
967 
968  void StdExpansion::v_NormVectorIProductWRTBase(const Array<OneD, const NekDouble> &Fx, const Array<OneD, const NekDouble> &Fy, Array< OneD, NekDouble> &outarray)
969  {
970  NEKERROR(ErrorUtil::efatal, "This function is not valid for this class");
971  }
972 
973  void StdExpansion::v_NormVectorIProductWRTBase(
974  const Array<OneD, const NekDouble> &Fx,
975  const Array<OneD, const NekDouble> &Fy,
976  const Array<OneD, const NekDouble> &Fz,
977  Array< OneD, NekDouble> &outarray)
978  {
979  NEKERROR(ErrorUtil::efatal, "This function is not valid for this class");
980  }
981 
982  DNekScalBlkMatSharedPtr StdExpansion::v_GetLocStaticCondMatrix(const LocalRegions::MatrixKey &mkey)
983  {
984  NEKERROR(ErrorUtil::efatal, "This function is only valid for LocalRegions");
985  return NullDNekScalBlkMatSharedPtr;
986  }
987 
988  void StdExpansion::v_DropLocStaticCondMatrix(const LocalRegions::MatrixKey &mkey)
989  {
990  NEKERROR(ErrorUtil::efatal, "This function is only valid for LocalRegions");
991  }
992 
993  StdRegions::Orientation StdExpansion::v_GetForient(int face)
994 
995  {
996  NEKERROR(ErrorUtil::efatal, "This function is only valid for three-dimensional LocalRegions");
997  return eDir1FwdDir1_Dir2FwdDir2;
998  }
999 
1000  StdRegions::Orientation StdExpansion::v_GetEorient(int edge)
1001  {
1002  NEKERROR(ErrorUtil::efatal, "This function is only valid for two-dimensional LocalRegions");
1003  return eForwards;
1004  }
1005 
1006  StdRegions::Orientation StdExpansion::v_GetPorient(int point)
1007  {
1008  NEKERROR(ErrorUtil::efatal, "This function is only valid for one-dimensional LocalRegions");
1009  return eFwd;
1010  }
1011 
1012 
1013  StdRegions::Orientation StdExpansion::v_GetCartesianEorient(int edge)
1014  {
1015  NEKERROR(ErrorUtil::efatal, "This function is only valid for two-dimensional LocalRegions");
1016  return eForwards;
1017  }
1018 
1019 
1020  void StdExpansion::v_SetCoeffsToOrientation(StdRegions::Orientation dir,
1021  Array<OneD, const NekDouble> &inarray,
1022  Array<OneD, NekDouble> &outarray)
1023  {
1024  NEKERROR(ErrorUtil::efatal, "This function is not defined for this shape");
1025  }
1026 
1027  void StdExpansion::v_SetCoeffsToOrientation(
1028  Array<OneD, NekDouble> &coeffs,
1029  StdRegions::Orientation dir)
1030  {
1031  NEKERROR(ErrorUtil::efatal, "This function is not defined for this shape");
1032  }
1033 
1034 
1035  NekDouble StdExpansion::v_StdPhysEvaluate(const Array<OneD, const NekDouble> &Lcoord,
1036  const Array<OneD, const NekDouble> &physvals)
1037 
1038  {
1039  NEKERROR(ErrorUtil::efatal, "This function is not defined for this shape");
1040  return 0;
1041  }
1042 
1043  void StdExpansion::v_LocCoordToLocCollapsed(const Array<OneD, const NekDouble>& xi,Array<OneD, NekDouble>& eta)
1044  {
1045  NEKERROR(ErrorUtil::efatal, "This function is not defined for this shape");
1046  }
1047 
1048  int StdExpansion::v_GetNedges() const
1049  {
1050  ASSERTL0(false, "This function is needs defining for this shape");
1051  return 0;
1052  }
1053 
1054  int StdExpansion::v_GetNfaces() const
1055  {
1056  ASSERTL0(false, "This function is needs defining for this shape");
1057  return 0;
1058  }
1059 
1060  int StdExpansion::v_NumBndryCoeffs() const
1061  {
1062  ASSERTL0(false, "This function is needs defining for this shape");
1063  return 0;
1064  }
1065 
1066  int StdExpansion::v_NumDGBndryCoeffs() const
1067  {
1068  ASSERTL0(false, "This function is needs defining for this shape");
1069  return 0;
1070  }
1071 
1072  int StdExpansion::v_GetEdgeNcoeffs(const int i) const
1073  {
1074  ASSERTL0(false, "This function is not valid or not defined");
1075  return 0;
1076  }
1077 
1078  int StdExpansion::v_GetTotalEdgeIntNcoeffs() const
1079  {
1080  ASSERTL0(false, "This function is not valid or not defined");
1081  return 0;
1082  }
1083 
1084  int StdExpansion::v_GetEdgeNumPoints(const int i) const
1085  {
1086  ASSERTL0(false, "This function is not valid or not defined");
1087  return 0;
1088  }
1089 
1090  int StdExpansion::v_DetCartesianDirOfEdge(const int edge)
1091  {
1092  ASSERTL0(false, "This function is not valid or not defined");
1093  return 0;
1094  }
1095 
1096  const LibUtilities::BasisKey StdExpansion::v_DetEdgeBasisKey(const int i) const
1097  {
1098  ASSERTL0(false, "This function is not valid or not defined");
1099  return LibUtilities::NullBasisKey;
1100  }
1101 
1102  const LibUtilities::BasisKey StdExpansion::v_DetFaceBasisKey(const int i, const int k) const
1103  {
1104  ASSERTL0(false, "This function is not valid or not defined");
1105  return LibUtilities::NullBasisKey;
1106  }
1107 
1108  int StdExpansion::v_GetFaceNumPoints(const int i) const
1109  {
1110  ASSERTL0(false, "This function is not valid or not defined");
1111  return 0;
1112  }
1113 
1114  int StdExpansion::v_GetFaceNcoeffs(const int i) const
1115  {
1116  ASSERTL0(false, "This function is not valid or not defined");
1117  return 0;
1118  }
1119 
1120  int StdExpansion::v_GetFaceIntNcoeffs(const int i) const
1121  {
1122  ASSERTL0(false, "This function is not valid or not defined");
1123  return 0;
1124  }
1125 
1126  int StdExpansion::v_GetTotalFaceIntNcoeffs() const
1127  {
1128  ASSERTL0(false, "This function is not valid or not defined");
1129  return 0;
1130  }
1131 
1132  LibUtilities::PointsKey StdExpansion::v_GetFacePointsKey(const int i, const int j) const
1133  {
1134  ASSERTL0(false, "This function is not valid or not defined");
1135  return LibUtilities::NullPointsKey;
1136  }
1137 
1138  LibUtilities::BasisType StdExpansion::v_GetEdgeBasisType(const int i) const
1139  {
1140  ASSERTL0(false, "This function is not valid or not defined");
1141 
1142  return LibUtilities::eNoBasisType;
1143  }
1144 
1145  LibUtilities::ShapeType StdExpansion::v_DetShapeType() const
1146  {
1147  ASSERTL0(false, "This expansion does not have a shape type defined");
1148  return LibUtilities::eNoShapeType;
1149  }
1150 
1151  int StdExpansion::v_GetShapeDimension() const
1152  {
1153  ASSERTL0(false, "This function is not valid or not defined");
1154  return 0;
1155  }
1156 
1157  bool StdExpansion::v_IsBoundaryInteriorExpansion()
1158  {
1159  ASSERTL0(false,"This function has not been defined for this expansion");
1160  return false;
1161  }
1162 
1163  void StdExpansion::v_IProductWRTDerivBase (const int dir,
1164  const Array<OneD, const NekDouble>& inarray,
1165  Array<OneD, NekDouble> &outarray)
1166  {
1167  NEKERROR(ErrorUtil::efatal, "This method has not been defined");
1168  }
1169 
1170  /**
1171  *
1172  */
1173  void StdExpansion::v_FwdTrans_BndConstrained(const Array<OneD, const NekDouble>& inarray,
1174  Array<OneD, NekDouble> &outarray)
1175  {
1176  NEKERROR(ErrorUtil::efatal, "This method has not been defined");
1177  }
1178 
1179 
1180  /**
1181  * @brief Integrates the specified function over the domain.
1182  * @see StdRegions#StdExpansion#Integral.
1183  */
1184  NekDouble StdExpansion::v_Integral(const Array<OneD, const NekDouble>& inarray )
1185  {
1186  NEKERROR(ErrorUtil::efatal, "This function is only valid for "
1187  "local expansions");
1188  return 0;
1189  }
1190 
1191 
1192  void StdExpansion::v_AddRobinMassMatrix(const int edgeid, const Array<OneD, const NekDouble > &primCoeffs, DNekMatSharedPtr &inoutmat)
1193  {
1194  NEKERROR(ErrorUtil::efatal, "This function is only valid for "
1195  "specific element types");
1196  }
1197 
1198  void StdExpansion::v_AddRobinEdgeContribution(const int edgeid, const Array<OneD, const NekDouble > &primCoeffs, Array<OneD, NekDouble> &coeffs)
1199  {
1200  NEKERROR(ErrorUtil::efatal, "This function is only valid for "
1201  "specific element types");
1202  }
1203 
1204  /**
1205  * @brief Calculate the derivative of the physical points
1206  * @see StdRegions#StdExpansion#PhysDeriv
1207  */
1208  void StdExpansion::v_PhysDeriv (const Array<OneD, const NekDouble>& inarray,
1209  Array<OneD, NekDouble> &out_d1,
1210  Array<OneD, NekDouble> &out_d2,
1211  Array<OneD, NekDouble> &out_d3)
1212  {
1213  NEKERROR(ErrorUtil::efatal, "This function is only valid for "
1214  "local expansions");
1215  }
1216 
1217  void StdExpansion::v_PhysDeriv_s(const Array<OneD, const NekDouble>& inarray,
1218  Array<OneD, NekDouble> &out_ds)
1219  {
1220  NEKERROR(ErrorUtil::efatal, "This function is only valid for "
1221  "local expansions");
1222  }
1223  void StdExpansion::v_PhysDeriv_n(const Array<OneD, const NekDouble>& inarray,
1224  Array<OneD, NekDouble>& out_dn)
1225  {
1226  NEKERROR(ErrorUtil::efatal, "This function is only valid for "
1227  "local expansions");
1228  }
1229 
1230  /**
1231  * @brief Calculate the derivative of the physical points in a
1232  * given direction
1233  * @see StdRegions#StdExpansion#PhysDeriv
1234  */
1235  void StdExpansion::v_PhysDeriv(const int dir,
1236  const Array<OneD, const NekDouble>& inarray,
1237  Array<OneD, NekDouble> &out_d0)
1238 
1239  {
1240  NEKERROR(ErrorUtil::efatal, "This function is only valid for "
1241  "specific element types");
1242  }
1243 
1244  /**
1245  * @brief Physical derivative along a direction vector.
1246  * @see StdRegions#StdExpansion#PhysDirectionalDeriv
1247  */
1248  void StdExpansion::v_PhysDirectionalDeriv(const Array<OneD, const NekDouble>& inarray,
1249  const Array<OneD, const NekDouble>& direction,
1250  Array<OneD, NekDouble> &outarray)
1251  {
1252  NEKERROR(ErrorUtil::efatal, "This function is only valid for "
1253  "specific element types");
1254  }
1255 
1256  void StdExpansion::v_StdPhysDeriv (const Array<OneD, const NekDouble>& inarray,
1257  Array<OneD, NekDouble> &out_d1,
1258  Array<OneD, NekDouble> &out_d2,
1259  Array<OneD, NekDouble> &out_d3)
1260  {
1261  NEKERROR(ErrorUtil::efatal, "Method does not exist for this shape");
1262  }
1263 
1264  void StdExpansion::v_StdPhysDeriv (const int dir,
1265  const Array<OneD, const NekDouble>& inarray,
1266  Array<OneD, NekDouble> &outarray)
1267  {
1268  NEKERROR(ErrorUtil::efatal, "Method does not exist for this shape");
1269  }
1270 
1271 
1272  NekDouble StdExpansion::v_PhysEvaluate(const Array<OneD, const NekDouble>& coords, const Array<OneD, const NekDouble>& physvals)
1273  {
1274  NEKERROR(ErrorUtil::efatal, "Method does not exist for this shape");
1275  return 0;
1276  }
1277 
1278 
1279  NekDouble StdExpansion::v_PhysEvaluate(const Array<OneD, DNekMatSharedPtr > & I, const Array<OneD, const NekDouble>& physvals)
1280  {
1281  NEKERROR(ErrorUtil::efatal, "Method does not exist for this shape");
1282  return 0;
1283  }
1284 
1285 
1286  void StdExpansion::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
1287  {
1288  NEKERROR(ErrorUtil::efatal, "This function has not "
1289  "been defined for this shape");
1290  }
1291 
1292  DNekMatSharedPtr StdExpansion::v_GenMatrix(const StdMatrixKey &mkey)
1293  {
1294  NEKERROR(ErrorUtil::efatal, "This function has not "
1295  "been defined for this element");
1296  DNekMatSharedPtr returnval;
1297  return returnval;
1298  }
1299 
1300  DNekMatSharedPtr StdExpansion::v_CreateStdMatrix(const StdMatrixKey &mkey)
1301  {
1302  NEKERROR(ErrorUtil::efatal, "This function has not "
1303  "been defined for this element");
1304  DNekMatSharedPtr returnval;
1305  return returnval;
1306  }
1307 
1308  void StdExpansion::v_GetCoords(Array<OneD, NekDouble> &coords_0,
1309  Array<OneD, NekDouble> &coords_1,
1310  Array<OneD, NekDouble> &coords_2)
1311  {
1312  NEKERROR(ErrorUtil::efatal, "Write coordinate definition method");
1313  }
1314 
1315  void StdExpansion::v_GetCoord(const Array<OneD, const NekDouble>& Lcoord,
1316  Array<OneD, NekDouble> &coord)
1317  {
1318  NEKERROR(ErrorUtil::efatal, "Write coordinate definition method");
1319  }
1320 
1321  int StdExpansion::v_GetCoordim(void)
1322  {
1323  NEKERROR(ErrorUtil::efatal, "Write method");
1324  return -1;
1325  }
1326 
1327  void StdExpansion::v_GetBoundaryMap(Array<OneD, unsigned int>& outarray)
1328  {
1329  NEKERROR(ErrorUtil::efatal,"Method does not exist for this shape" );
1330  }
1331 
1332  void StdExpansion::v_GetInteriorMap(Array<OneD, unsigned int>& outarray)
1333  {
1334  NEKERROR(ErrorUtil::efatal,"Method does not exist for this shape" );
1335  }
1336 
1337  int StdExpansion::v_GetVertexMap(const int localVertexId,
1338  bool useCoeffPacking)
1339  {
1340  NEKERROR(ErrorUtil::efatal,"Method does not exist for this shape" );
1341  return 0;
1342  }
1343 
1344  void StdExpansion::v_GetEdgeInteriorMap(const int eid, const Orientation edgeOrient,
1345  Array<OneD, unsigned int> &maparray,
1346  Array<OneD, int> &signarray)
1347  {
1348  NEKERROR(ErrorUtil::efatal,"Method does not exist for this shape" );
1349  }
1350 
1351  void StdExpansion::v_GetFaceInteriorMap(const int fid, const Orientation faceOrient,
1352  Array<OneD, unsigned int> &maparray,
1353  Array<OneD, int> &signarray)
1354  {
1355  NEKERROR(ErrorUtil::efatal,"Method does not exist for this shape" );
1356  }
1357 
1358  void StdExpansion::v_GetEdgeToElementMap(const int eid, const Orientation edgeOrient,
1359  Array<OneD, unsigned int> &maparray,
1360  Array<OneD, int> &signarray)
1361  {
1362  NEKERROR(ErrorUtil::efatal,"Method does not exist for this shape" );
1363  }
1364 
1365  void StdExpansion::v_GetFaceToElementMap(const int fid, const Orientation faceOrient,
1366  Array<OneD, unsigned int> &maparray,
1367  Array<OneD, int> &signarray,
1368  int nummodesA, int nummodesB)
1369  {
1370  NEKERROR(ErrorUtil::efatal,"Method does not exist for this shape" );
1371  }
1372 
1373  void StdExpansion::v_GetEdgePhysVals(const int edge, const Array<OneD, const NekDouble> &inarray, Array<OneD,NekDouble> &outarray)
1374  {
1375  NEKERROR(ErrorUtil::efatal,"Method does not exist for this shape or library" );
1376  }
1377 
1378  void StdExpansion::v_GetEdgePhysVals(const int edge, const boost::shared_ptr<StdExpansion> &EdgeExp, const Array<OneD, const NekDouble> &inarray, Array<OneD,NekDouble> &outarray)
1379  {
1380  NEKERROR(ErrorUtil::efatal,"Method does not exist for this shape or library" );
1381  }
1382 
1383  void StdExpansion::v_GetTracePhysVals(const int edge, const boost::shared_ptr<StdExpansion> &EdgeExp, const Array<OneD, const NekDouble> &inarray, Array<OneD,NekDouble> &outarray, StdRegions::Orientation orient)
1384  {
1385  NEKERROR(ErrorUtil::efatal,"Method does not exist for this shape or library" );
1386  }
1387 
1388  void StdExpansion::v_GetVertexPhysVals(const int vertex, const Array<OneD, const NekDouble> &inarray, NekDouble &outarray)
1389  {
1390  NEKERROR(ErrorUtil::efatal,"Method does not exist for this shape or library" );
1391  }
1392 
1393  void StdExpansion::v_GetEdgeInterpVals(const int edge,const Array<OneD, const NekDouble> &inarray,Array<OneD,NekDouble> &outarray)
1394  {
1395  NEKERROR(ErrorUtil::efatal,"Method does not exist for this shape or library" );
1396  }
1397 
1398  void StdExpansion::v_GetEdgeQFactors(
1399  const int edge,
1400  Array<OneD, NekDouble> &outarray)
1401  {
1402  NEKERROR(ErrorUtil::efatal,
1403  "Method does not exist for this shape or library");
1404  }
1405 
1406  void StdExpansion::v_GetFacePhysVals(
1407  const int face,
1408  const boost::shared_ptr<StdExpansion> &FaceExp,
1409  const Array<OneD, const NekDouble> &inarray,
1410  Array<OneD, NekDouble> &outarray,
1411  StdRegions::Orientation orient)
1412  {
1413  NEKERROR(ErrorUtil::efatal,"Method does not exist for this shape or library" );
1414  }
1415 
1416  void StdExpansion::v_MultiplyByQuadratureMetric(
1417  const Array<OneD, const NekDouble> &inarray,
1418  Array<OneD, NekDouble> &outarray)
1419  {
1420  v_MultiplyByStdQuadratureMetric(inarray, outarray);
1421  }
1422 
1423  void StdExpansion::v_MultiplyByStdQuadratureMetric(
1424  const Array<OneD, const NekDouble> &inarray,
1425  Array<OneD, NekDouble> &outarray)
1426  {
1427  NEKERROR(ErrorUtil::efatal, "Method does not exist for this shape or library");
1428  }
1429 
1430  const boost::shared_ptr<SpatialDomains::GeomFactors>& StdExpansion::v_GetMetricInfo() const
1431  {
1432  NEKERROR(ErrorUtil::efatal, "This function is only valid for LocalRegions");
1433  return SpatialDomains::NullGeomFactorsSharedPtr;
1434 
1435  }
1436 
1437  void StdExpansion::v_BwdTrans_SumFac(const Array<OneD, const NekDouble>& inarray,
1438  Array<OneD, NekDouble> &outarray)
1439  {
1440  NEKERROR(ErrorUtil::efatal,"Method does not exist for this shape" );
1441  }
1442 
1443  void StdExpansion::v_IProductWRTBase_SumFac(const Array<OneD, const NekDouble>& inarray,
1444  Array<OneD, NekDouble> &outarray)
1445  {
1446  NEKERROR(ErrorUtil::efatal,"Method does not exist for this shape" );
1447  }
1448 
1449  void StdExpansion::v_IProductWRTDerivBase_SumFac(const int dir,
1450  const Array<OneD, const NekDouble>& inarray,
1451  Array<OneD, NekDouble> &outarray)
1452  {
1453  NEKERROR(ErrorUtil::efatal,"Method does not exist for this shape" );
1454  }
1455 
1456  void StdExpansion::v_MassMatrixOp(const Array<OneD, const NekDouble> &inarray,
1457  Array<OneD,NekDouble> &outarray,
1458  const StdMatrixKey &mkey)
1459  {
1460  // If this function is not reimplemented on shape level, the function
1461  // below will be called
1462  MassMatrixOp_MatFree(inarray,outarray,mkey);
1463  }
1464 
1465  void StdExpansion::v_LaplacianMatrixOp(const Array<OneD, const NekDouble> &inarray,
1466  Array<OneD,NekDouble> &outarray,
1467  const StdMatrixKey &mkey)
1468  {
1469  // If this function is not reimplemented on shape level, the function
1470  // below will be called
1471  LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
1472  }
1473 
1474 
1475  void StdExpansion::v_SVVLaplacianFilter(Array<OneD,NekDouble> &array,
1476  const StdMatrixKey &mkey)
1477  {
1478  ASSERTL0(false, "This function is not defined in StdExpansion.");
1479  }
1480 
1481  void StdExpansion::v_ReduceOrderCoeffs(int numMin,
1482  const Array<OneD, const NekDouble> &inarray,
1483  Array<OneD, NekDouble> &outarray)
1484  {
1485  ASSERTL0(false, "This function is not defined in StdExpansion.");
1486  }
1487 
1488  void StdExpansion::v_LaplacianMatrixOp(const int k1, const int k2,
1489  const Array<OneD, const NekDouble> &inarray,
1490  Array<OneD,NekDouble> &outarray,
1491  const StdMatrixKey &mkey)
1492  {
1493  // If this function is not reimplemented on shape level, the function
1494  // below will be called
1495  LaplacianMatrixOp_MatFree(k1,k2,inarray,outarray,mkey);
1496  }
1497 
1498  void StdExpansion::v_WeakDerivMatrixOp(const int i,
1499  const Array<OneD, const NekDouble> &inarray,
1500  Array<OneD,NekDouble> &outarray,
1501  const StdMatrixKey &mkey)
1502  {
1503  // If this function is not reimplemented on shape level, the function
1504  // below will be called
1505  WeakDerivMatrixOp_MatFree(i,inarray,outarray,mkey);
1506 
1507  }
1508 
1509  void StdExpansion::v_WeakDirectionalDerivMatrixOp(const Array<OneD, const NekDouble> &inarray,
1510  Array<OneD,NekDouble> &outarray,
1511  const StdMatrixKey &mkey)
1512  {
1513  // If this function is not reimplemented on shape level, the function
1514  // below will be called
1515  WeakDirectionalDerivMatrixOp_MatFree(inarray,outarray,mkey);
1516 
1517  }
1518 
1519  void StdExpansion::v_MassLevelCurvatureMatrixOp(const Array<OneD, const NekDouble> &inarray,
1520  Array<OneD,NekDouble> &outarray,
1521  const StdMatrixKey &mkey)
1522  {
1523  // If this function is not reimplemented on shape level, the function
1524  // below will be called
1525  MassLevelCurvatureMatrixOp_MatFree(inarray,outarray,mkey);
1526  }
1527 
1528  void StdExpansion::v_LinearAdvectionDiffusionReactionMatrixOp(const Array<OneD,
1529  const NekDouble> &inarray,
1530  Array<OneD,NekDouble> &outarray,
1531  const StdMatrixKey &mkey, bool addDiffusionTerm)
1532  {
1533  // If this function is not reimplemented on shape level, the function
1534  // below will be called
1535  LinearAdvectionDiffusionReactionMatrixOp_MatFree(inarray,outarray,mkey,addDiffusionTerm);
1536 
1537  }
1538 
1539  void StdExpansion::v_HelmholtzMatrixOp(const Array<OneD, const NekDouble> &inarray,
1540  Array<OneD,NekDouble> &outarray,
1541  const StdMatrixKey &mkey)
1542  {
1543  // If this function is not reimplemented on shape level, the function
1544  // below will be called
1545  HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
1546  }
1547 
1548  void StdExpansion::v_LaplacianMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
1549  Array<OneD,NekDouble> &outarray,
1550  const StdMatrixKey &mkey)
1551  {
1552  // If this function is not reimplemented on shape level, the function
1553  // below will be called
1554  LaplacianMatrixOp_MatFree_GenericImpl(inarray,outarray,mkey);
1555  }
1556 
1557  void StdExpansion::v_LaplacianMatrixOp_MatFree_Kernel(
1558  const Array<OneD, const NekDouble> &inarray,
1559  Array<OneD, NekDouble> &outarray,
1560  Array<OneD, NekDouble> &wsp)
1561  {
1562  ASSERTL0(false, "Not implemented.");
1563  }
1564 
1565  void StdExpansion::v_HelmholtzMatrixOp_MatFree(const Array<OneD, const NekDouble> &inarray,
1566  Array<OneD,NekDouble> &outarray,
1567  const StdMatrixKey &mkey)
1568  {
1569  // If this function is not reimplemented on shape level, the function
1570  // below will be called
1571  HelmholtzMatrixOp_MatFree_GenericImpl(inarray,outarray,mkey);
1572  }
1573 
1574  const NormalVector & StdExpansion::v_GetEdgeNormal(const int edge) const
1575  {
1576  ASSERTL0(false, "Cannot get edge normals for this expansion.");
1577  static NormalVector result;
1578  return result;
1579  }
1580 
1581  void StdExpansion::v_ComputeEdgeNormal(const int edge)
1582  {
1583  ASSERTL0(false, "Cannot compute edge normal for this expansion.");
1584  }
1585 
1586  void StdExpansion::v_NegateEdgeNormal(const int edge)
1587  {
1588  ASSERTL0(false, "Not implemented.");
1589  }
1590 
1591  bool StdExpansion::v_EdgeNormalNegated(const int edge)
1592  {
1593  ASSERTL0(false, "Not implemented.");
1594  return false;
1595  }
1596 
1597  void StdExpansion::v_ComputeFaceNormal(const int face)
1598  {
1599  ASSERTL0(false, "Cannot compute face normal for this expansion.");
1600  }
1601 
1602  void StdExpansion::v_NegateFaceNormal(const int face)
1603  {
1604  ASSERTL0(false, "Not implemented.");
1605  }
1606 
1607  void StdExpansion::v_ComputeVertexNormal(const int vertex)
1608  {
1609  ASSERTL0(false, "Cannot compute vertex normal for this expansion.");
1610  }
1611 
1612  const NormalVector & StdExpansion::v_GetFaceNormal(const int face) const
1613  {
1614  ASSERTL0(false, "Cannot get face normals for this expansion.");
1615  static NormalVector result;
1616  return result;
1617  }
1618 
1619  const NormalVector & StdExpansion::v_GetVertexNormal(const int vertex) const
1620  {
1621  ASSERTL0(false, "Cannot get vertex normals for this expansion.");
1622  static NormalVector result;
1623  return result;
1624  }
1625 
1626  const NormalVector & StdExpansion::v_GetSurfaceNormal(const int id) const
1627  {
1628  ASSERTL0(false, "Cannot get face normals for this expansion.");
1629  static NormalVector result;
1630  return result;
1631  }
1632 
1633  Array<OneD, unsigned int>
1634  StdExpansion::v_GetEdgeInverseBoundaryMap(int eid)
1635  {
1636  ASSERTL0(false, "Not implemented.");
1637  Array<OneD, unsigned int> noinversemap(1);
1638  return noinversemap;
1639  }
1640 
1641  Array<OneD, unsigned int>
1642  StdExpansion::v_GetFaceInverseBoundaryMap(int fid,
1643  StdRegions::Orientation faceOrient)
1644  {
1645  ASSERTL0(false, "Not implemented.");
1646  Array<OneD, unsigned int> noinversemap(1);
1647  return noinversemap;
1648  }
1649 
1650  DNekMatSharedPtr
1651  StdExpansion::v_BuildInverseTransformationMatrix(
1652  const DNekScalMatSharedPtr & m_transformationmatrix)
1653  {
1654  NEKERROR(ErrorUtil::efatal, "This function is only valid for LocalRegions");
1655  return NullDNekMatSharedPtr;
1656  }
1657 
1658  void StdExpansion::PhysInterpToSimplexEquiSpaced(
1659  const Array<OneD, const NekDouble> &inarray,
1660  Array<OneD, NekDouble> &outarray)
1661  {
1662  LibUtilities::ShapeType shape = DetShapeType();
1663  StdMatrixKey Ikey(ePhysInterpToEquiSpaced, shape, *this);
1664  DNekMatSharedPtr intmat = GetStdMatrix(Ikey);
1665 
1666  int nqtot = 1;
1667  int nqbase;
1668  int np = 0;
1669  for(int i = 0; i < m_base.num_elements(); ++i)
1670  {
1671  nqbase = m_base[i]->GetNumPoints();
1672  nqtot *= nqbase;
1673  np = max(np,nqbase);
1674  }
1675 
1676  NekVector<NekDouble> in (nqtot,inarray,eWrapper);
1677  NekVector<NekDouble> out(LibUtilities::GetNumberOfCoefficients(shape,np,np,np),outarray,eWrapper);
1678  out = (*intmat) * in;
1679  }
1680 
1681  void StdExpansion::v_GetSimplexEquiSpacedConnectivity(
1682  Array<OneD, int> &conn,
1683  bool standard)
1684  {
1685  ASSERTL0(false, "Not implemented.");
1686  }
1687  }//end namespace
1688 }//end namespace
Generated on Wed Mar 4 2015 16:33:17 for Nektar++ by   doxygen 1.8.1.2