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