Nektar++
TriExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File TriExp.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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 // The above copyright notice and this permission notice shall be included
20 // in all copies or substantial portions of the Software.
21 //
22 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 // DEALINGS IN THE SOFTWARE.
29 //
30 // Description: Expasion for triangular elements.
31 //
32 ///////////////////////////////////////////////////////////////////////////////
33 
34 #include <boost/core/ignore_unused.hpp>
35 
36 #include <LibUtilities/Foundations/InterpCoeff.h>
37 #include <LocalRegions/TriExp.h>
38 #include <LocalRegions/SegExp.h>
39 #include <LocalRegions/Expansion3D.h>
40 #include <StdRegions/StdNodalTriExp.h>
41 #include <LibUtilities/Foundations/Interp.h>
42 
43 using namespace std;
44 
45 namespace Nektar
46 {
47  namespace LocalRegions
48  {
49  TriExp::TriExp(const LibUtilities::BasisKey &Ba,
50  const LibUtilities::BasisKey &Bb,
51  const SpatialDomains::TriGeomSharedPtr &geom):
52  StdExpansion (LibUtilities::StdTriData::getNumberOfCoefficients(Ba.GetNumModes(),(Bb.GetNumModes())),2,Ba,Bb),
53  StdExpansion2D(LibUtilities::StdTriData::getNumberOfCoefficients(Ba.GetNumModes(),(Bb.GetNumModes())),Ba,Bb),
54  StdTriExp(Ba,Bb),
55  Expansion (geom),
56  Expansion2D (geom),
57  m_matrixManager(
58  std::bind(&TriExp::CreateMatrix, this, std::placeholders::_1),
59  std::string("TriExpMatrix")),
60  m_staticCondMatrixManager(
61  std::bind(&TriExp::CreateStaticCondMatrix, this, std::placeholders::_1),
62  std::string("TriExpStaticCondMatrix"))
63  {
64  }
65 
66 
67  TriExp::TriExp(const TriExp &T):
68  StdExpansion(T),
69  StdExpansion2D(T),
70  StdTriExp(T),
71  Expansion(T),
72  Expansion2D(T),
73  m_matrixManager(T.m_matrixManager),
74  m_staticCondMatrixManager(T.m_staticCondMatrixManager)
75  {
76  }
77 
78 
79  TriExp::~TriExp()
80  {
81  }
82 
83 
84 
85  NekDouble TriExp::v_Integral(const Array<OneD, const NekDouble> &inarray)
86  {
87  int nquad0 = m_base[0]->GetNumPoints();
88  int nquad1 = m_base[1]->GetNumPoints();
89  Array<OneD, const NekDouble> jac = m_metricinfo->GetJac(GetPointsKeys());
90  NekDouble ival;
91  Array<OneD,NekDouble> tmp(nquad0*nquad1);
92 
93  // multiply inarray with Jacobian
94  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
95  {
96  Vmath::Vmul(nquad0*nquad1, jac, 1, inarray, 1,tmp, 1);
97  }
98  else
99  {
100  Vmath::Smul(nquad0*nquad1, jac[0], inarray, 1, tmp, 1);
101  }
102 
103  // call StdQuadExp version;
104  ival = StdTriExp::v_Integral(tmp);
105  return ival;
106  }
107 
108  void TriExp::v_PhysDeriv(const Array<OneD, const NekDouble> & inarray,
109  Array<OneD,NekDouble> &out_d0,
110  Array<OneD,NekDouble> &out_d1,
111  Array<OneD,NekDouble> &out_d2)
112  {
113  int nquad0 = m_base[0]->GetNumPoints();
114  int nquad1 = m_base[1]->GetNumPoints();
115  int nqtot = nquad0*nquad1;
116  const Array<TwoD, const NekDouble>& df
117  = m_metricinfo->GetDerivFactors(GetPointsKeys());
118 
119  Array<OneD,NekDouble> diff0(2*nqtot);
120  Array<OneD,NekDouble> diff1(diff0+nqtot);
121 
122  StdTriExp::v_PhysDeriv(inarray, diff0, diff1);
123 
124  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
125  {
126  if(out_d0.num_elements())
127  {
128  Vmath::Vmul (nqtot,df[0],1,diff0,1, out_d0, 1);
129  Vmath::Vvtvp (nqtot,df[1],1,diff1,1, out_d0, 1, out_d0,1);
130  }
131 
132  if(out_d1.num_elements())
133  {
134  Vmath::Vmul (nqtot,df[2],1,diff0,1, out_d1, 1);
135  Vmath::Vvtvp (nqtot,df[3],1,diff1,1, out_d1, 1, out_d1,1);
136  }
137 
138  if(out_d2.num_elements())
139  {
140  Vmath::Vmul (nqtot,df[4],1,diff0,1, out_d2, 1);
141  Vmath::Vvtvp (nqtot,df[5],1,diff1,1, out_d2, 1, out_d2,1);
142  }
143  }
144  else // regular geometry
145  {
146  if(out_d0.num_elements())
147  {
148  Vmath::Smul (nqtot, df[0][0], diff0, 1, out_d0, 1);
149  Blas::Daxpy (nqtot, df[1][0], diff1, 1, out_d0, 1);
150  }
151 
152  if(out_d1.num_elements())
153  {
154  Vmath::Smul (nqtot, df[2][0], diff0, 1, out_d1, 1);
155  Blas::Daxpy (nqtot, df[3][0], diff1, 1, out_d1, 1);
156  }
157 
158  if(out_d2.num_elements())
159  {
160  Vmath::Smul (nqtot, df[4][0], diff0, 1, out_d2, 1);
161  Blas::Daxpy (nqtot, df[5][0], diff1, 1, out_d2, 1);
162  }
163  }
164  }
165 
166 
167  void TriExp::v_PhysDeriv(const int dir,
168  const Array<OneD, const NekDouble>& inarray,
169  Array<OneD, NekDouble> &outarray)
170  {
171  switch(dir)
172  {
173  case 0:
174  {
175  PhysDeriv(inarray, outarray, NullNekDouble1DArray, NullNekDouble1DArray);
176  }
177  break;
178  case 1:
179  {
180  PhysDeriv(inarray, NullNekDouble1DArray, outarray, NullNekDouble1DArray);
181  }
182  break;
183  case 2:
184  {
185  PhysDeriv(inarray, NullNekDouble1DArray, NullNekDouble1DArray, outarray);
186  }
187  break;
188  default:
189  {
190  ASSERTL1(false,"input dir is out of range");
191  }
192  break;
193  }
194  }
195 
196  void TriExp::v_PhysDirectionalDeriv(
197  const Array<OneD, const NekDouble> &inarray,
198  const Array<OneD, const NekDouble> &direction,
199  Array<OneD, NekDouble> &out)
200  {
201  if(! out.num_elements())
202  {
203  return;
204  }
205 
206  int nquad0 = m_base[0]->GetNumPoints();
207  int nquad1 = m_base[1]->GetNumPoints();
208  int nqtot = nquad0*nquad1;
209 
210  const Array<TwoD, const NekDouble>& df =
211  m_metricinfo->GetDerivFactors(GetPointsKeys());
212 
213  Array<OneD,NekDouble> diff0(2*nqtot);
214  Array<OneD,NekDouble> diff1(diff0+nqtot);
215 
216  // diff0 = du/d_xi, diff1 = du/d_eta
217  StdTriExp::v_PhysDeriv(inarray, diff0, diff1);
218 
219  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
220  {
221  Array<OneD, Array<OneD, NekDouble> > tangmat(2);
222 
223 
224  // D^v_xi = v_x*d_xi/dx + v_y*d_xi/dy + v_z*d_xi/dz
225  // D^v_eta = v_x*d_eta/dx + v_y*d_eta/dy + v_z*d_eta/dz
226  for (int i=0; i< 2; ++i)
227  {
228  tangmat[i] = Array<OneD, NekDouble>(nqtot,0.0);
229  for (int k=0; k<(m_geom->GetCoordim()); ++k)
230  {
231  Vmath::Vvtvp(nqtot,&df[2*k+i][0],1,&direction[k*nqtot],1,&tangmat[i][0],1,&tangmat[i][0],1);
232  }
233  }
234 
235  /// D_v = D^v_xi * du/d_xi + D^v_eta * du/d_eta
236  Vmath::Vmul (nqtot,&tangmat[0][0],1,&diff0[0],1, &out[0], 1);
237  Vmath::Vvtvp (nqtot,&tangmat[1][0],1,&diff1[0],1, &out[0], 1, &out[0],1);
238  }
239  else
240  {
241  Array<OneD, Array<OneD, NekDouble> > tangmat(2);
242 
243  for (int i=0; i< 2; ++i)
244  {
245  tangmat[i] = Array<OneD, NekDouble>(nqtot,0.0);
246  for (int k=0; k<(m_geom->GetCoordim()); ++k)
247  {
248  Vmath::Svtvp(nqtot, df[2*k+i][0],
249  &direction[k*nqtot], 1,
250  &tangmat[i][0], 1, &tangmat[i][0], 1);
251  }
252  }
253 
254  /// D_v = D^v_xi * du/d_xi + D^v_eta * du/d_eta
255  Vmath::Vmul (nqtot, &tangmat[0][0], 1,
256  &diff0[0], 1,
257  &out[0], 1);
258 
259  Vmath::Vvtvp (nqtot, &tangmat[1][0], 1,
260  &diff1[0], 1,
261  &out[0], 1,
262  &out[0], 1);
263  }
264  }
265 
266 
267  void TriExp::v_FwdTrans(const Array<OneD, const NekDouble> & inarray,
268  Array<OneD,NekDouble> &outarray)
269  {
270  IProductWRTBase(inarray,outarray);
271 
272  // get Mass matrix inverse
273  MatrixKey masskey(StdRegions::eInvMass,
274  DetShapeType(),*this);
275  DNekScalMatSharedPtr matsys = m_matrixManager[masskey];
276 
277  // copy inarray in case inarray == outarray
278  NekVector<NekDouble> in (m_ncoeffs,outarray,eCopy);
279  NekVector<NekDouble> out(m_ncoeffs,outarray,eWrapper);
280 
281  out = (*matsys)*in;
282  }
283 
284 
285  void TriExp::v_FwdTrans_BndConstrained(const Array<OneD, const NekDouble>& inarray,
286  Array<OneD, NekDouble> &outarray)
287  {
288  int i,j;
289  int npoints[2] = {m_base[0]->GetNumPoints(),
290  m_base[1]->GetNumPoints()};
291  int nmodes[2] = {m_base[0]->GetNumModes(),
292  m_base[1]->GetNumModes()};
293 
294  fill(outarray.get(), outarray.get()+m_ncoeffs, 0.0 );
295 
296  if(nmodes[0] == 1 && nmodes[1] == 1)
297  {
298  outarray[0] = inarray[0];
299  return;
300  }
301 
302  Array<OneD, NekDouble> physEdge[3];
303  Array<OneD, NekDouble> coeffEdge[3];
304  for(i = 0; i < 3; i++)
305  {
306  // define physEdge and add 1 so can interpolate grl10 points if necessary
307  physEdge[i] = Array<OneD, NekDouble>(max(npoints[i!=0],npoints[0]));
308  coeffEdge[i] = Array<OneD, NekDouble>(nmodes[i!=0]);
309  }
310 
311  for(i = 0; i < npoints[0]; i++)
312  {
313  physEdge[0][i] = inarray[i];
314  }
315 
316  // extract data in cartesian directions
317  for(i = 0; i < npoints[1]; i++)
318  {
319  physEdge[1][i] = inarray[npoints[0]-1+i*npoints[0]];
320  physEdge[2][i] = inarray[i*npoints[0]];
321  }
322 
323  SegExpSharedPtr segexp[3];
324  segexp[0] = MemoryManager<LocalRegions::SegExp>::AllocateSharedPtr(m_base[0]->GetBasisKey(),GetGeom2D()->GetEdge(0));
325 
326  if(m_base[1]->GetPointsType() == LibUtilities::eGaussLobattoLegendre)
327  {
328  for(i = 1; i < 3; i++)
329  {
330  segexp[i] = MemoryManager<LocalRegions::SegExp>::AllocateSharedPtr(m_base[i!=0]->GetBasisKey(),GetGeom2D()->GetEdge(i));
331  }
332  }
333  else // interploate using edge 0 GLL distribution
334  {
335  for(i = 1; i < 3; i++)
336  {
337  segexp[i] = MemoryManager<LocalRegions::SegExp>::AllocateSharedPtr(m_base[0]->GetBasisKey(),GetGeom2D()->GetEdge(i));
338 
339  LibUtilities::Interp1D(m_base[1]->GetPointsKey(),physEdge[i],
340  m_base[0]->GetPointsKey(),physEdge[i]);
341  }
342  npoints[1] = npoints[0];
343  }
344 
345 
346  Array<OneD, unsigned int> mapArray;
347  Array<OneD, int> signArray;
348  NekDouble sign;
349  // define an orientation to get EdgeToElmtMapping from Cartesian data
350  StdRegions::Orientation orient[3] = {StdRegions::eForwards,StdRegions::eForwards,
351  StdRegions::eForwards};
352 
353  for(i = 0; i < 3; i++)
354  {
355  segexp[i]->FwdTrans_BndConstrained(physEdge[i],coeffEdge[i]);
356 
357  // this orient goes with the one above and so could
358  // probably set both to eForwards
359  GetEdgeToElementMap(i,orient[i],mapArray,signArray);
360  for(j=0; j < nmodes[i!=0]; j++)
361  {
362  sign = (NekDouble) signArray[j];
363  outarray[ mapArray[j] ] = sign * coeffEdge[i][j];
364  }
365  }
366 
367  int nBoundaryDofs = NumBndryCoeffs();
368  int nInteriorDofs = m_ncoeffs - nBoundaryDofs;
369 
370  if (nInteriorDofs > 0) {
371  Array<OneD, NekDouble> tmp0(m_ncoeffs);
372  Array<OneD, NekDouble> tmp1(m_ncoeffs);
373 
374  StdRegions::StdMatrixKey stdmasskey(StdRegions::eMass,DetShapeType(),*this);
375  MassMatrixOp(outarray,tmp0,stdmasskey);
376  IProductWRTBase(inarray,tmp1);
377 
378  Vmath::Vsub(m_ncoeffs, tmp1, 1, tmp0, 1, tmp1, 1);
379 
380  // get Mass matrix inverse (only of interior DOF)
381  // use block (1,1) of the static condensed system
382  // note: this block alreay contains the inverse matrix
383  MatrixKey masskey(StdRegions::eMass,DetShapeType(),*this);
384  DNekScalMatSharedPtr matsys = (m_staticCondMatrixManager[masskey])->GetBlock(1,1);
385 
386  Array<OneD, NekDouble> rhs(nInteriorDofs);
387  Array<OneD, NekDouble> result(nInteriorDofs);
388 
389  GetInteriorMap(mapArray);
390 
391  for(i = 0; i < nInteriorDofs; i++)
392  {
393  rhs[i] = tmp1[ mapArray[i] ];
394  }
395 
396  Blas::Dgemv('N', nInteriorDofs, nInteriorDofs, matsys->Scale(), &((matsys->GetOwnedMatrix())->GetPtr())[0],
397  nInteriorDofs,rhs.get(),1,0.0,result.get(),1);
398 
399  for(i = 0; i < nInteriorDofs; i++)
400  {
401  outarray[ mapArray[i] ] = result[i];
402  }
403  }
404  }
405 
406 
407  void TriExp::v_IProductWRTBase(const Array<OneD, const NekDouble>& inarray,
408  Array<OneD, NekDouble> &outarray)
409  {
410  IProductWRTBase_SumFac(inarray,outarray);
411  }
412 
413 
414  void TriExp::v_IProductWRTDerivBase(const int dir,
415  const Array<OneD, const NekDouble>& inarray,
416  Array<OneD, NekDouble> & outarray)
417  {
418  IProductWRTDerivBase_SumFac(dir,inarray,outarray);
419  }
420 
421 
422  void TriExp::v_IProductWRTBase_SumFac(const Array<OneD, const NekDouble>& inarray,
423  Array<OneD, NekDouble> &outarray,
424  bool multiplybyweights)
425  {
426  int nquad0 = m_base[0]->GetNumPoints();
427  int nquad1 = m_base[1]->GetNumPoints();
428  int order0 = m_base[0]->GetNumModes();
429 
430  if(multiplybyweights)
431  {
432  Array<OneD,NekDouble> tmp(nquad0*nquad1+nquad1*order0);
433  Array<OneD,NekDouble> wsp(tmp+nquad0*nquad1);
434 
435  MultiplyByQuadratureMetric(inarray,tmp);
436  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),m_base[1]->GetBdata(),tmp,outarray,wsp);
437  }
438  else
439  {
440  Array<OneD,NekDouble> wsp(+nquad1*order0);
441 
442  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),m_base[1]->GetBdata(),
443  inarray,outarray,wsp);
444  }
445  }
446 
447 
448  void TriExp::v_IProductWRTBase_MatOp(const Array<OneD, const NekDouble>& inarray,
449  Array<OneD, NekDouble> &outarray)
450  {
451  int nq = GetTotPoints();
452  MatrixKey iprodmatkey(StdRegions::eIProductWRTBase,DetShapeType(),*this);
453  DNekScalMatSharedPtr iprodmat = m_matrixManager[iprodmatkey];
454 
455  Blas::Dgemv('N',m_ncoeffs,nq,iprodmat->Scale(),(iprodmat->GetOwnedMatrix())->GetPtr().get(),
456  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
457 
458  }
459 
460 
461  void TriExp::v_IProductWRTDerivBase_SumFac(const int dir,
462  const Array<OneD, const NekDouble>& inarray,
463  Array<OneD, NekDouble> & outarray)
464  {
465  ASSERTL1((dir==0)||(dir==1)||(dir==2),"Invalid direction.");
466  ASSERTL1((dir==2)?(m_geom->GetCoordim()==3):true,"Invalid direction.");
467 
468  int i;
469  int nquad0 = m_base[0]->GetNumPoints();
470  int nquad1 = m_base[1]->GetNumPoints();
471  int nqtot = nquad0*nquad1;
472  int nmodes0 = m_base[0]->GetNumModes();
473  int wspsize = max(max(nqtot,m_ncoeffs),nquad1*nmodes0);
474 
475  const Array<TwoD, const NekDouble>& df =
476  m_metricinfo->GetDerivFactors(GetPointsKeys());
477 
478  Array<OneD, NekDouble> tmp0 (6*wspsize);
479  Array<OneD, NekDouble> tmp1 (tmp0 + wspsize);
480  Array<OneD, NekDouble> tmp2 (tmp0 + 2*wspsize);
481  Array<OneD, NekDouble> tmp3 (tmp0 + 3*wspsize);
482  Array<OneD, NekDouble> gfac0(tmp0 + 4*wspsize);
483  Array<OneD, NekDouble> gfac1(tmp0 + 5*wspsize);
484 
485  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
486  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
487 
488  // set up geometric factor: 2/(1-z1)
489  for(i = 0; i < nquad1; ++i)
490  {
491  gfac0[i] = 2.0/(1-z1[i]);
492  }
493  for(i = 0; i < nquad0; ++i)
494  {
495  gfac1[i] = 0.5*(1+z0[i]);
496  }
497 
498  for(i = 0; i < nquad1; ++i)
499  {
500  Vmath::Smul(nquad0,gfac0[i],&inarray[0]+i*nquad0,1,&tmp0[0]+i*nquad0,1);
501  }
502 
503  for(i = 0; i < nquad1; ++i)
504  {
505  Vmath::Vmul(nquad0,&gfac1[0],1,&tmp0[0]+i*nquad0,1,&tmp1[0]+i*nquad0,1);
506  }
507 
508  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
509  {
510  Vmath::Vmul(nqtot,&df[2*dir][0], 1,&tmp0[0], 1,&tmp0[0],1);
511  Vmath::Vmul(nqtot,&df[2*dir+1][0],1,&tmp1[0], 1,&tmp1[0],1);
512  Vmath::Vmul(nqtot,&df[2*dir+1][0],1,&inarray[0],1,&tmp2[0],1);
513  }
514  else
515  {
516  Vmath::Smul(nqtot, df[2*dir][0], tmp0, 1, tmp0, 1);
517  Vmath::Smul(nqtot, df[2*dir+1][0], tmp1, 1, tmp1, 1);
518  Vmath::Smul(nqtot, df[2*dir+1][0], inarray, 1, tmp2, 1);
519  }
520  Vmath::Vadd(nqtot, tmp0, 1, tmp1, 1, tmp1, 1);
521 
522  MultiplyByQuadratureMetric(tmp1,tmp1);
523  MultiplyByQuadratureMetric(tmp2,tmp2);
524 
525  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),m_base[1]->GetBdata() ,tmp1,tmp3 ,tmp0);
526  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata() ,m_base[1]->GetDbdata(),tmp2,outarray,tmp0);
527  Vmath::Vadd(m_ncoeffs, tmp3, 1, outarray, 1, outarray, 1);
528  }
529 
530 
531  void TriExp::v_IProductWRTDerivBase_MatOp(const int dir,
532  const Array<OneD, const NekDouble>& inarray,
533  Array<OneD, NekDouble> &outarray)
534  {
535  int nq = GetTotPoints();
536  StdRegions::MatrixType mtype = StdRegions::eIProductWRTDerivBase0;
537 
538  switch(dir)
539  {
540  case 0:
541  {
542  mtype = StdRegions::eIProductWRTDerivBase0;
543  }
544  break;
545  case 1:
546  {
547  mtype = StdRegions::eIProductWRTDerivBase1;
548  }
549  break;
550  case 2:
551  {
552  mtype = StdRegions::eIProductWRTDerivBase2;
553  }
554  break;
555  default:
556  {
557  ASSERTL1(false,"input dir is out of range");
558  }
559  break;
560  }
561 
562  MatrixKey iprodmatkey(mtype,DetShapeType(),*this);
563  DNekScalMatSharedPtr iprodmat = m_matrixManager[iprodmatkey];
564 
565  Blas::Dgemv('N',m_ncoeffs,nq,iprodmat->Scale(),(iprodmat->GetOwnedMatrix())->GetPtr().get(),
566  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
567 
568  }
569 
570 
571  void TriExp::v_IProductWRTDirectionalDerivBase(
572  const Array<OneD, const NekDouble>& direction,
573  const Array<OneD, const NekDouble>& inarray,
574  Array<OneD, NekDouble> & outarray)
575  {
576  IProductWRTDirectionalDerivBase_SumFac(direction,
577  inarray,outarray);
578  }
579 
580 
581  /**
582  * @brief Directinoal Derivative in the modal space in the dir
583  * direction of varcoeffs.
584  */
585  void TriExp::v_IProductWRTDirectionalDerivBase_SumFac(
586  const Array<OneD, const NekDouble>& direction,
587  const Array<OneD, const NekDouble>& inarray,
588  Array<OneD, NekDouble> & outarray)
589  {
590  int i;
591  int shapedim = 2;
592  int nquad0 = m_base[0]->GetNumPoints();
593  int nquad1 = m_base[1]->GetNumPoints();
594  int nqtot = nquad0*nquad1;
595  int nmodes0 = m_base[0]->GetNumModes();
596  int wspsize = max(max(nqtot,m_ncoeffs),nquad1*nmodes0);
597 
598  const Array<TwoD, const NekDouble>& df =
599  m_metricinfo->GetDerivFactors(GetPointsKeys());
600 
601  Array<OneD, NekDouble> tmp0 (6*wspsize);
602  Array<OneD, NekDouble> tmp1 (tmp0 + wspsize);
603  Array<OneD, NekDouble> tmp2 (tmp0 + 2*wspsize);
604  Array<OneD, NekDouble> tmp3 (tmp0 + 3*wspsize);
605  Array<OneD, NekDouble> gfac0(tmp0 + 4*wspsize);
606  Array<OneD, NekDouble> gfac1(tmp0 + 5*wspsize);
607 
608  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
609  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
610 
611  // set up geometric factor: 2/(1-z1)
612  for(i = 0; i < nquad1; ++i)
613  {
614  gfac0[i] = 2.0/(1-z1[i]);
615  }
616  for(i = 0; i < nquad0; ++i)
617  {
618  gfac1[i] = 0.5*(1+z0[i]);
619  }
620  for(i = 0; i < nquad1; ++i)
621  {
622  Vmath::Smul(nquad0, gfac0[i], &inarray[0] + i*nquad0, 1,
623  &tmp0[0] + i*nquad0, 1);
624  }
625  for(i = 0; i < nquad1; ++i)
626  {
627  Vmath::Vmul(nquad0, &gfac1[0], 1,
628  &tmp0[0] + i*nquad0, 1,
629  &tmp1[0] + i*nquad0, 1);
630  }
631 
632  // Compute gmat \cdot e^j
633  Array<OneD, Array<OneD, NekDouble> > dfdir(shapedim);
634  Expansion::ComputeGmatcdotMF(df, direction, dfdir);
635 
636  Vmath::Vmul(nqtot, &dfdir[0][0], 1, &tmp0[0], 1, &tmp0[0], 1);
637  Vmath::Vmul(nqtot, &dfdir[1][0], 1, &tmp1[0], 1, &tmp1[0], 1);
638  Vmath::Vmul(nqtot, &dfdir[1][0], 1, &inarray[0], 1, &tmp2[0], 1);
639 
640  Vmath::Vadd(nqtot, &tmp0[0], 1, &tmp1[0], 1, &tmp1[0], 1);
641 
642  MultiplyByQuadratureMetric(tmp1,tmp1);
643  MultiplyByQuadratureMetric(tmp2,tmp2);
644 
645  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
646  m_base[1]->GetBdata(),
647  tmp1, tmp3, tmp0);
648  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
649  m_base[1]->GetDbdata(),
650  tmp2, outarray, tmp0);
651  Vmath::Vadd(m_ncoeffs, tmp3, 1, outarray, 1, outarray, 1);
652  }
653 
654 
655  void TriExp::v_NormVectorIProductWRTBase(
656  const Array<OneD, const NekDouble> &Fx,
657  const Array<OneD, const NekDouble> &Fy,
658  const Array<OneD, const NekDouble> &Fz,
659  Array<OneD, NekDouble> &outarray)
660  {
661  int nq = m_base[0]->GetNumPoints()*m_base[1]->GetNumPoints();
662  Array<OneD, NekDouble > Fn(nq);
663 
664  const Array<OneD, const Array<OneD, NekDouble> > &normals =
665  GetLeftAdjacentElementExp()->GetFaceNormal(
666  GetLeftAdjacentElementFace());
667 
668  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
669  {
670  Vmath::Vvtvvtp(nq,&normals[0][0],1,&Fx[0],1,
671  &normals[1][0],1,&Fy[0],1,&Fn[0],1);
672  Vmath::Vvtvp (nq,&normals[2][0],1,&Fz[0],1,&Fn[0],1,&Fn[0],1);
673  }
674  else
675  {
676  Vmath::Svtsvtp(nq,normals[0][0],&Fx[0],1,
677  normals[1][0],&Fy[0],1,&Fn[0],1);
678  Vmath::Svtvp (nq,normals[2][0],&Fz[0],1,&Fn[0],1,&Fn[0],1);
679  }
680 
681  IProductWRTBase(Fn,outarray);
682  }
683 
684  void TriExp::v_NormVectorIProductWRTBase(
685  const Array<OneD, const Array<OneD, NekDouble> > &Fvec,
686  Array<OneD, NekDouble> &outarray)
687  {
688  NormVectorIProductWRTBase(Fvec[0], Fvec[1], Fvec[2], outarray);
689  }
690 
691  StdRegions::StdExpansionSharedPtr TriExp::v_GetStdExp(void) const
692  {
693 
694  return MemoryManager<StdRegions::StdTriExp>
695  ::AllocateSharedPtr(m_base[0]->GetBasisKey(),
696  m_base[1]->GetBasisKey());
697  }
698 
699  StdRegions::StdExpansionSharedPtr TriExp::v_GetLinStdExp(void) const
700  {
701  LibUtilities::BasisKey bkey0(m_base[0]->GetBasisType(),
702  2, m_base[0]->GetPointsKey());
703  LibUtilities::BasisKey bkey1(m_base[1]->GetBasisType(),
704  2, m_base[1]->GetPointsKey());
705 
706  return MemoryManager<StdRegions::StdTriExp>
707  ::AllocateSharedPtr( bkey0, bkey1);
708  }
709 
710  void TriExp::v_GetCoord(const Array<OneD, const NekDouble> &Lcoords,
711  Array<OneD,NekDouble> &coords)
712  {
713  int i;
714 
715  ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[1] <= 1.0 &&
716  Lcoords[1] >= -1.0 && Lcoords[1] <=1.0,
717  "Local coordinates are not in region [-1,1]");
718 
719  m_geom->FillGeom();
720 
721  for(i = 0; i < m_geom->GetCoordim(); ++i)
722  {
723  coords[i] = m_geom->GetCoord(i,Lcoords);
724  }
725  }
726 
727  void TriExp::v_GetCoords(
728  Array<OneD, NekDouble> &coords_0,
729  Array<OneD, NekDouble> &coords_1,
730  Array<OneD, NekDouble> &coords_2)
731  {
732  Expansion::v_GetCoords(coords_0, coords_1, coords_2);
733  }
734 
735 
736  /**
737  * Given the local cartesian coordinate \a Lcoord evaluate the
738  * value of physvals at this point by calling through to the
739  * StdExpansion method
740  */
741  NekDouble TriExp::v_StdPhysEvaluate(
742  const Array<OneD, const NekDouble> &Lcoord,
743  const Array<OneD, const NekDouble> &physvals)
744  {
745  // Evaluate point in local (eta) coordinates.
746  return StdTriExp::v_PhysEvaluate(Lcoord,physvals);
747  }
748 
749  NekDouble TriExp::v_PhysEvaluate(const Array<OneD, const NekDouble> &coord, const Array<OneD, const NekDouble> & physvals)
750  {
751  Array<OneD,NekDouble> Lcoord = Array<OneD,NekDouble>(2);
752 
753  ASSERTL0(m_geom,"m_geom not defined");
754  m_geom->GetLocCoords(coord,Lcoord);
755 
756  return StdTriExp::v_PhysEvaluate(Lcoord, physvals);
757  }
758 
759 
760  void TriExp::v_GetTracePhysVals(
761  const int edge,
762  const StdRegions::StdExpansionSharedPtr &EdgeExp,
763  const Array<OneD, const NekDouble> &inarray,
764  Array<OneD,NekDouble> &outarray,
765  StdRegions::Orientation orient)
766  {
767  boost::ignore_unused(orient);
768  v_GetEdgePhysVals(edge,EdgeExp,inarray,outarray);
769  }
770 
771  void TriExp::v_GetEdgePhysVals(
772  const int edge,
773  const Array<OneD, const NekDouble> &inarray,
774  Array<OneD,NekDouble> &outarray)
775  {
776  int nquad0 = m_base[0]->GetNumPoints();
777  int nquad1 = m_base[1]->GetNumPoints();
778 
779  StdRegions::Orientation edgedir = GetEorient(edge);
780  switch(edge)
781  {
782  case 0:
783  if (edgedir == StdRegions::eForwards)
784  {
785  Vmath::Vcopy(nquad0,&(inarray[0]),1,&(outarray[0]),1);
786  }
787  else
788  {
789  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0-1),-1,
790  &(outarray[0]),1);
791  }
792  break;
793  case 1:
794  if (edgedir == StdRegions::eForwards)
795  {
796  Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0-1),nquad0,
797  &(outarray[0]),1);
798  }
799  else
800  {
801  Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0*nquad1-1),
802  -nquad0, &(outarray[0]),1);
803  }
804  break;
805  case 2:
806  if (edgedir == StdRegions::eForwards)
807  {
808  Vmath::Vcopy(nquad1,&(inarray[0]) + nquad0*(nquad1-1),
809  -nquad0,&(outarray[0]),1);
810  }
811  else
812  {
813  Vmath::Vcopy(nquad1,&(inarray[0]),nquad0,
814  &(outarray[0]),1);
815  }
816  break;
817  default:
818  ASSERTL0(false,"edge value (< 3) is out of range");
819  break;
820  }
821  }
822 
823  void TriExp::v_GetEdgePhysVals(
824  const int edge,
825  const StdRegions::StdExpansionSharedPtr &EdgeExp,
826  const Array<OneD, const NekDouble> &inarray,
827  Array<OneD, NekDouble> &outarray)
828  {
829  int nquad0 = m_base[0]->GetNumPoints();
830  int nquad1 = m_base[1]->GetNumPoints();
831 
832  // Extract in Cartesian direction because we have to deal with
833  // e.g. Gauss-Radau points.
834  switch(edge)
835  {
836  case 0:
837  Vmath::Vcopy(nquad0, &(inarray[0]), 1, &(outarray[0]), 1);
838  break;
839  case 1:
840  Vmath::Vcopy(nquad1, &(inarray[0])+(nquad0-1),
841  nquad0, &(outarray[0]), 1);
842  break;
843  case 2:
844  Vmath::Vcopy(nquad1, &(inarray[0]), nquad0, &(outarray[0]), 1);
845  break;
846  default:
847  ASSERTL0(false,"edge value (< 3) is out of range");
848  break;
849  }
850 
851  ASSERTL1(EdgeExp->GetBasis(0)->GetPointsType()
852  == LibUtilities::eGaussLobattoLegendre,
853  "Edge expansion should be GLL");
854 
855  // Interpolate if required
856  if(m_base[edge?1:0]->GetPointsKey() != EdgeExp->GetBasis(0)->GetPointsKey())
857  {
858  Array<OneD,NekDouble> outtmp(max(nquad0,nquad1));
859 
860  outtmp = outarray;
861 
862  LibUtilities::Interp1D(m_base[edge?1:0]->GetPointsKey(),
863  outtmp,
864  EdgeExp->GetBasis(0)->GetPointsKey(),
865  outarray);
866  }
867 
868  StdRegions::Orientation orient = GetEorient(edge);
869 
870  //Reverse data if necessary
871  if(orient == StdRegions::eBackwards)
872  {
873  Vmath::Reverse(EdgeExp->GetNumPoints(0),&outarray[0],1,
874  &outarray[0],1);
875  }
876 
877  }
878 
879 
880  void TriExp::v_GetEdgeInterpVals(
881  const int edge,const Array<OneD, const NekDouble> &inarray,
882  Array<OneD, NekDouble> &outarray)
883  {
884  boost::ignore_unused(edge, inarray, outarray);
885  ASSERTL0(false,
886  "Routine not implemented for triangular elements");
887  }
888 
889  void TriExp::v_GetEdgeQFactors(
890  const int edge,
891  Array<OneD, NekDouble> &outarray)
892  {
893  boost::ignore_unused(edge, outarray);
894  ASSERTL0(false,
895  "Routine not implemented for triangular elements");
896  }
897 
898 
899 
900  void TriExp::v_GetEdgePhysMap(
901  const int edge,
902  Array<OneD, int> &outarray)
903  {
904  int nquad0 = m_base[0]->GetNumPoints();
905  int nquad1 = m_base[1]->GetNumPoints();
906 
907  // Get points in Cartesian orientation
908  switch (edge)
909  {
910  case 0:
911  outarray = Array<OneD, int>(nquad0);
912  for (int i = 0; i < nquad0; ++i)
913  {
914  outarray[i] = i;
915  }
916  break;
917  case 1:
918  outarray = Array<OneD, int>(nquad1);
919  for (int i = 0; i < nquad1; ++i)
920  {
921  outarray[i] = (nquad0-1) + i * nquad0;
922  }
923  break;
924  case 2:
925  outarray = Array<OneD, int>(nquad1);
926  for (int i = 0; i < nquad1; ++i)
927  {
928  outarray[i] = i*nquad0;
929  }
930  break;
931  default:
932  ASSERTL0(false, "edge value (< 3) is out of range");
933  break;
934  }
935 
936  }
937 
938 
939  void TriExp::v_ComputeEdgeNormal(const int edge)
940  {
941  int i;
942  const SpatialDomains::GeomFactorsSharedPtr & geomFactors = GetGeom()->GetMetricInfo();
943 
944  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
945  for(i = 0; i < ptsKeys.size(); ++i)
946  {
947  // Need at least 2 points for computing normals
948  if (ptsKeys[i].GetNumPoints() == 1)
949  {
950  LibUtilities::PointsKey pKey(2, ptsKeys[i].GetPointsType());
951  ptsKeys[i] = pKey;
952  }
953  }
954 
955  const SpatialDomains::GeomType type = geomFactors->GetGtype();
956  const Array<TwoD, const NekDouble> & df = geomFactors->GetDerivFactors(ptsKeys);
957  const Array<OneD, const NekDouble> & jac = geomFactors->GetJac(ptsKeys);
958  int nqe = m_base[0]->GetNumPoints();
959  int dim = GetCoordim();
960 
961  m_edgeNormals[edge] = Array<OneD, Array<OneD, NekDouble> >(dim);
962  Array<OneD, Array<OneD, NekDouble> > &normal = m_edgeNormals[edge];
963  for (i = 0; i < dim; ++i)
964  {
965  normal[i] = Array<OneD, NekDouble>(nqe);
966  }
967 
968  // Regular geometry case
969  if((type == SpatialDomains::eRegular)||(type == SpatialDomains::eMovingRegular))
970  {
971  NekDouble fac;
972  // Set up normals
973  switch(edge)
974  {
975  case 0:
976  for(i = 0; i < GetCoordim(); ++i)
977  {
978  Vmath::Fill(nqe,-df[2*i+1][0],normal[i],1);
979  }
980  break;
981  case 1:
982  for(i = 0; i < GetCoordim(); ++i)
983  {
984  Vmath::Fill(nqe,df[2*i+1][0] + df[2*i][0],normal[i],1);
985  }
986  break;
987  case 2:
988  for(i = 0; i < GetCoordim(); ++i)
989  {
990  Vmath::Fill(nqe,-df[2*i][0],normal[i],1);
991  }
992  break;
993  default:
994  ASSERTL0(false,"Edge is out of range (edge < 3)");
995  }
996 
997  // normalise
998  fac = 0.0;
999  for(i =0 ; i < GetCoordim(); ++i)
1000  {
1001  fac += normal[i][0]*normal[i][0];
1002  }
1003  fac = 1.0/sqrt(fac);
1004  for (i = 0; i < GetCoordim(); ++i)
1005  {
1006  Vmath::Smul(nqe,fac,normal[i],1,normal[i],1);
1007  }
1008  }
1009  else // Set up deformed normals
1010  {
1011  int j;
1012 
1013  int nquad0 = ptsKeys[0].GetNumPoints();
1014  int nquad1 = ptsKeys[1].GetNumPoints();
1015 
1016  LibUtilities::PointsKey from_key;
1017 
1018  Array<OneD,NekDouble> normals(GetCoordim()*max(nquad0,nquad1),0.0);
1019  Array<OneD,NekDouble> edgejac(GetCoordim()*max(nquad0,nquad1),0.0);
1020 
1021  // Extract Jacobian along edges and recover local
1022  // derivates (dx/dr) for polynomial interpolation by
1023  // multiplying m_gmat by jacobian
1024  switch(edge)
1025  {
1026  case 0:
1027  for(j = 0; j < nquad0; ++j)
1028  {
1029  edgejac[j] = jac[j];
1030  for(i = 0; i < GetCoordim(); ++i)
1031  {
1032  normals[i*nquad0+j] = -df[2*i+1][j]*edgejac[j];
1033  }
1034  }
1035  from_key = ptsKeys[0];
1036  break;
1037  case 1:
1038  for(j = 0; j < nquad1; ++j)
1039  {
1040  edgejac[j] = jac[nquad0*j+nquad0-1];
1041  for(i = 0; i < GetCoordim(); ++i)
1042  {
1043  normals[i*nquad1+j] = (df[2*i][nquad0*j + nquad0-1] + df[2*i+1][nquad0*j + nquad0-1])*edgejac[j];
1044  }
1045  }
1046  from_key = ptsKeys[1];
1047  break;
1048  case 2:
1049  for(j = 0; j < nquad1; ++j)
1050  {
1051  edgejac[j] = jac[nquad0*j];
1052  for(i = 0; i < GetCoordim(); ++i)
1053  {
1054  normals[i*nquad1+j] = -df[2*i][nquad0*j]*edgejac[j];
1055  }
1056  }
1057  from_key = ptsKeys[1];
1058  break;
1059  default:
1060  ASSERTL0(false,"edge is out of range (edge < 3)");
1061 
1062  }
1063 
1064  int nq = from_key.GetNumPoints();
1065  Array<OneD,NekDouble> work(nqe,0.0);
1066 
1067  // interpolate Jacobian and invert
1068  LibUtilities::Interp1D(from_key,jac,m_base[0]->GetPointsKey(),work);
1069  Vmath::Sdiv(nqe,1.0,&work[0],1,&work[0],1);
1070 
1071  // interpolate
1072  for(i = 0; i < GetCoordim(); ++i)
1073  {
1074  LibUtilities::Interp1D(from_key,&normals[i*nq],m_base[0]->GetPointsKey(),&normal[i][0]);
1075  Vmath::Vmul(nqe,work,1,normal[i],1,normal[i],1);
1076  }
1077 
1078  //normalise normal vectors
1079  Vmath::Zero(nqe,work,1);
1080  for(i = 0; i < GetCoordim(); ++i)
1081  {
1082  Vmath::Vvtvp(nqe,normal[i],1, normal[i],1,work,1,work,1);
1083  }
1084 
1085  Vmath::Vsqrt(nqe,work,1,work,1);
1086  Vmath::Sdiv(nqe,1.0,work,1,work,1);
1087 
1088  for(i = 0; i < GetCoordim(); ++i)
1089  {
1090  Vmath::Vmul(nqe,normal[i],1,work,1,normal[i],1);
1091  }
1092  }
1093  if(GetGeom()->GetEorient(edge) == StdRegions::eBackwards)
1094  {
1095  for(i = 0; i < GetCoordim(); ++i)
1096  {
1097  if(geomFactors->GetGtype() == SpatialDomains::eDeformed)
1098  {
1099  Vmath::Reverse(nqe, normal[i], 1, normal[i],1);
1100  }
1101  }
1102  }
1103  }
1104 
1105  int TriExp::v_GetCoordim()
1106  {
1107  return m_geom->GetCoordim();
1108  }
1109 
1110 
1111  void TriExp::v_ExtractDataToCoeffs(
1112  const NekDouble *data,
1113  const std::vector<unsigned int > &nummodes,
1114  const int mode_offset,
1115  NekDouble * coeffs,
1116  std::vector<LibUtilities::BasisType> &fromType)
1117  {
1118  boost::ignore_unused(fromType);
1119 
1120  int data_order0 = nummodes[mode_offset];
1121  int fillorder0 = min(m_base[0]->GetNumModes(),data_order0);
1122  int data_order1 = nummodes[mode_offset+1];
1123  int order1 = m_base[1]->GetNumModes();
1124  int fillorder1 = min(order1,data_order1);
1125 
1126  switch(m_base[0]->GetBasisType())
1127  {
1128  case LibUtilities::eModified_A:
1129  {
1130  int i;
1131  int cnt = 0;
1132  int cnt1 = 0;
1133 
1134  ASSERTL1(m_base[1]->GetBasisType() == LibUtilities::eModified_B,
1135  "Extraction routine not set up for this basis");
1136 
1137  Vmath::Zero(m_ncoeffs,coeffs,1);
1138  for(i = 0; i < fillorder0; ++i)
1139  {
1140  Vmath::Vcopy(fillorder1-i,&data[cnt],1,&coeffs[cnt1],1);
1141  cnt += data_order1-i;
1142  cnt1 += order1-i;
1143  }
1144  }
1145  break;
1146  default:
1147  ASSERTL0(false,"basis is either not set up or not hierarchicial");
1148  }
1149  }
1150 
1151 
1152  StdRegions::Orientation TriExp::v_GetEorient(int edge)
1153  {
1154  return GetGeom2D()->GetEorient(edge);
1155  }
1156 
1157  const LibUtilities::BasisSharedPtr& TriExp::v_GetBasis(int dir) const
1158  {
1159  ASSERTL1(dir >= 0 &&dir <= 1,"input dir is out of range");
1160  return m_base[dir];
1161  }
1162 
1163 
1164  int TriExp::v_GetNumPoints(const int dir) const
1165  {
1166  return GetNumPoints(dir);
1167  }
1168 
1169 
1170  DNekMatSharedPtr TriExp::v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
1171  {
1172  DNekMatSharedPtr returnval;
1173  switch(mkey.GetMatrixType())
1174  {
1175  case StdRegions::eHybridDGHelmholtz:
1176  case StdRegions::eHybridDGLamToU:
1177  case StdRegions::eHybridDGLamToQ0:
1178  case StdRegions::eHybridDGLamToQ1:
1179  case StdRegions::eHybridDGLamToQ2:
1180  case StdRegions::eHybridDGHelmBndLam:
1181  case StdRegions::eInvLaplacianWithUnityMean:
1182  returnval = Expansion2D::v_GenMatrix(mkey);
1183  break;
1184  default:
1185  returnval = StdTriExp::v_GenMatrix(mkey);
1186  break;
1187  }
1188 
1189  return returnval;
1190  }
1191 
1192 
1193  DNekMatSharedPtr TriExp::v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
1194  {
1195  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
1196  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
1197  StdRegions::StdTriExpSharedPtr tmp = MemoryManager<StdTriExp>::
1198  AllocateSharedPtr(bkey0,bkey1);
1199 
1200  return tmp->GetStdMatrix(mkey);
1201  }
1202 
1203 
1204  DNekScalMatSharedPtr TriExp::CreateMatrix(const MatrixKey &mkey)
1205  {
1206  DNekScalMatSharedPtr returnval;
1207  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
1208 
1209  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
1210 
1211  switch(mkey.GetMatrixType())
1212  {
1213  case StdRegions::eMass:
1214  {
1215  if((m_metricinfo->GetGtype() == SpatialDomains::eDeformed)||
1216  (mkey.GetNVarCoeff()))
1217  {
1218  NekDouble one = 1.0;
1219  DNekMatSharedPtr mat = GenMatrix(mkey);
1220  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1221  }
1222  else
1223  {
1224  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1225  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1226  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
1227  }
1228  }
1229  break;
1230  case StdRegions::eInvMass:
1231  {
1232  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1233  {
1234  NekDouble one = 1.0;
1235  StdRegions::StdMatrixKey masskey(StdRegions::eMass,DetShapeType(),
1236  *this);
1237  DNekMatSharedPtr mat = GenMatrix(masskey);
1238  mat->Invert();
1239 
1240  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1241  }
1242  else
1243  {
1244  NekDouble fac = 1.0/(m_metricinfo->GetJac(ptsKeys))[0];
1245  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1246  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(fac,mat);
1247 
1248  }
1249  }
1250  break;
1251  case StdRegions::eWeakDeriv0:
1252  case StdRegions::eWeakDeriv1:
1253  case StdRegions::eWeakDeriv2:
1254  {
1255  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed || mkey.GetNVarCoeff())
1256  {
1257  NekDouble one = 1.0;
1258  DNekMatSharedPtr mat = GenMatrix(mkey);
1259 
1260  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1261  }
1262  else
1263  {
1264  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1265  Array<TwoD, const NekDouble> df = m_metricinfo->GetDerivFactors(ptsKeys);
1266  int dir = 0;
1267  switch(mkey.GetMatrixType())
1268  {
1269  case StdRegions::eWeakDeriv0:
1270  dir = 0;
1271  break;
1272  case StdRegions::eWeakDeriv1:
1273  dir = 1;
1274  break;
1275  case StdRegions::eWeakDeriv2:
1276  dir = 2;
1277  break;
1278  default:
1279  break;
1280  }
1281 
1282  MatrixKey deriv0key(StdRegions::eWeakDeriv0,
1283  mkey.GetShapeType(), *this);
1284  MatrixKey deriv1key(StdRegions::eWeakDeriv1,
1285  mkey.GetShapeType(), *this);
1286 
1287  DNekMat &deriv0 = *GetStdMatrix(deriv0key);
1288  DNekMat &deriv1 = *GetStdMatrix(deriv1key);
1289 
1290  int rows = deriv0.GetRows();
1291  int cols = deriv1.GetColumns();
1292 
1293  DNekMatSharedPtr WeakDeriv = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols);
1294  (*WeakDeriv) = df[2*dir][0]*deriv0 + df[2*dir+1][0]*deriv1;
1295 
1296  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,WeakDeriv);
1297  }
1298  }
1299  break;
1300  case StdRegions::eWeakDirectionalDeriv:
1301  {
1302  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1303  mkey.GetNVarCoeff())
1304  {
1305  NekDouble one = 1.0;
1306  DNekMatSharedPtr mat = GenMatrix(mkey);
1307 
1308  returnval = MemoryManager<DNekScalMat>::
1309  AllocateSharedPtr(one,mat);
1310  }
1311  else
1312  {
1313  int shapedim = 2;
1314 
1315  // dfdirxi = tan_{xi_x} * d \xi/dx
1316  // + tan_{xi_y} * d \xi/dy
1317  // + tan_{xi_z} * d \xi/dz
1318  // dfdireta = tan_{eta_x} * d \eta/dx
1319  // + tan_{xi_y} * d \xi/dy
1320  // + tan_{xi_z} * d \xi/dz
1321  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1322  Array<TwoD, const NekDouble> df =
1323  m_metricinfo->GetDerivFactors(ptsKeys);
1324 
1325  Array<OneD, NekDouble> direction =
1326  mkey.GetVarCoeff(StdRegions::eVarCoeffMF);
1327 
1328  // d / dx = df[0]*deriv0 + df[1]*deriv1
1329  // d / dy = df[2]*deriv0 + df[3]*deriv1
1330  // d / dz = df[4]*deriv0 + df[5]*deriv1
1331 
1332  // dfdir[dir] = e \cdot (d/dx, d/dy, d/dz)
1333  // = (e^0 * df[0] + e^1 * df[2]
1334  // + e^2 * df[4]) * deriv0
1335  // + (e^0 * df[1] + e^1 * df[3]
1336  // + e^2 * df[5]) * deriv1
1337  // dfdir[dir] = e^0 * df[2 * 0 + dir]
1338  // + e^1 * df[2 * 1 + dir]
1339  // + e^2 * df [ 2 * 2 + dir]
1340  Array<OneD, Array<OneD, NekDouble> > dfdir(shapedim);
1341  Expansion::ComputeGmatcdotMF(df, direction, dfdir);
1342 
1343  StdRegions::VarCoeffMap dfdirxi;
1344  StdRegions::VarCoeffMap dfdireta;
1345 
1346  dfdirxi[StdRegions::eVarCoeffWeakDeriv] = dfdir[0];
1347  dfdireta[StdRegions::eVarCoeffWeakDeriv] = dfdir[1];
1348 
1349  MatrixKey derivxikey( StdRegions::eWeakDeriv0,
1350  mkey.GetShapeType(), *this,
1351  StdRegions::NullConstFactorMap,
1352  dfdirxi);
1353  MatrixKey derivetakey( StdRegions::eWeakDeriv1,
1354  mkey.GetShapeType(), *this,
1355  StdRegions::NullConstFactorMap,
1356  dfdireta);
1357 
1358  DNekMat &derivxi = *GetStdMatrix(derivxikey);
1359  DNekMat &deriveta = *GetStdMatrix(derivetakey);
1360 
1361  int rows = derivxi.GetRows();
1362  int cols = deriveta.GetColumns();
1363 
1364  DNekMatSharedPtr WeakDirDeriv = MemoryManager<DNekMat>::
1365  AllocateSharedPtr(rows,cols);
1366 
1367  (*WeakDirDeriv) = derivxi + deriveta;
1368 
1369  // Add (\nabla \cdot e^k ) Mass
1370  StdRegions::VarCoeffMap DiveMass;
1371  DiveMass[StdRegions::eVarCoeffMass] =
1372  mkey.GetVarCoeff(StdRegions::eVarCoeffMFDiv);
1373  StdRegions::StdMatrixKey stdmasskey(
1374  StdRegions::eMass,
1375  mkey.GetShapeType(), *this,
1376  StdRegions::NullConstFactorMap,
1377  DiveMass);
1378 
1379  DNekMatSharedPtr DiveMassmat = GetStdMatrix(stdmasskey);
1380 
1381  (*WeakDirDeriv) = (*WeakDirDeriv) + (*DiveMassmat);
1382 
1383  returnval = MemoryManager<DNekScalMat>::
1384  AllocateSharedPtr(jac,WeakDirDeriv);
1385  }
1386  break;
1387  }
1388  case StdRegions::eLaplacian:
1389  {
1390  if( (m_metricinfo->GetGtype() == SpatialDomains::eDeformed) ||
1391  (mkey.GetNVarCoeff() > 0)||(mkey.ConstFactorExists(StdRegions::eFactorSVVCutoffRatio)))
1392  {
1393  NekDouble one = 1.0;
1394  DNekMatSharedPtr mat = GenMatrix(mkey);
1395 
1396  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1397  }
1398  else
1399  {
1400  MatrixKey lap00key(StdRegions::eLaplacian00,
1401  mkey.GetShapeType(), *this);
1402  MatrixKey lap01key(StdRegions::eLaplacian01,
1403  mkey.GetShapeType(), *this);
1404  MatrixKey lap11key(StdRegions::eLaplacian11,
1405  mkey.GetShapeType(), *this);
1406 
1407  DNekMat &lap00 = *GetStdMatrix(lap00key);
1408  DNekMat &lap01 = *GetStdMatrix(lap01key);
1409  DNekMat &lap11 = *GetStdMatrix(lap11key);
1410 
1411  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1412  Array<TwoD, const NekDouble> gmat =
1413  m_metricinfo->GetGmat(ptsKeys);
1414 
1415  int rows = lap00.GetRows();
1416  int cols = lap00.GetColumns();
1417 
1418  DNekMatSharedPtr lap = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols);
1419 
1420  (*lap) = gmat[0][0] * lap00 +
1421  gmat[1][0] * (lap01 + Transpose(lap01)) +
1422  gmat[3][0] * lap11;
1423 
1424  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,lap);
1425  }
1426  }
1427  break;
1428  case StdRegions::eInvLaplacianWithUnityMean:
1429  {
1430  DNekMatSharedPtr mat = GenMatrix(mkey);
1431  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0,mat);
1432  }
1433  break;
1434  case StdRegions::eHelmholtz:
1435  {
1436  NekDouble factor = mkey.GetConstFactor(StdRegions::eFactorLambda);
1437 
1438  MatrixKey masskey(mkey, StdRegions::eMass);
1439  DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
1440 
1441  MatrixKey lapkey(mkey, StdRegions::eLaplacian);
1442  DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
1443 
1444  int rows = LapMat.GetRows();
1445  int cols = LapMat.GetColumns();
1446 
1447  DNekMatSharedPtr helm = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols);
1448 
1449  NekDouble one = 1.0;
1450  (*helm) = LapMat + factor*MassMat;
1451 
1452  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,helm);
1453  }
1454  break;
1455  case StdRegions::eIProductWRTBase:
1456  {
1457  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1458  {
1459  NekDouble one = 1.0;
1460  DNekMatSharedPtr mat = GenMatrix(mkey);
1461  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1462  }
1463  else
1464  {
1465  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1466  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1467  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
1468  }
1469  }
1470  break;
1471  case StdRegions::eIProductWRTDerivBase0:
1472  case StdRegions::eIProductWRTDerivBase1:
1473  case StdRegions::eIProductWRTDerivBase2:
1474  {
1475  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1476  {
1477  NekDouble one = 1.0;
1478  DNekMatSharedPtr mat = GenMatrix(mkey);
1479  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1480  }
1481  else
1482  {
1483  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1484 
1485  const Array<TwoD, const NekDouble>& df = m_metricinfo->GetDerivFactors(ptsKeys);
1486  int dir = 0;
1487 
1488  switch(mkey.GetMatrixType())
1489  {
1490  case StdRegions::eIProductWRTDerivBase0:
1491  dir = 0;
1492  break;
1493  case StdRegions::eIProductWRTDerivBase1:
1494  dir = 1;
1495  break;
1496  case StdRegions::eIProductWRTDerivBase2:
1497  dir = 2;
1498  break;
1499  default:
1500  break;
1501  }
1502 
1503  MatrixKey iProdDeriv0Key(StdRegions::eIProductWRTDerivBase0,
1504  mkey.GetShapeType(), *this);
1505  MatrixKey iProdDeriv1Key(StdRegions::eIProductWRTDerivBase1,
1506  mkey.GetShapeType(), *this);
1507 
1508  DNekMat &stdiprod0 = *GetStdMatrix(iProdDeriv0Key);
1509  DNekMat &stdiprod1 = *GetStdMatrix(iProdDeriv0Key);
1510 
1511  int rows = stdiprod0.GetRows();
1512  int cols = stdiprod1.GetColumns();
1513 
1514  DNekMatSharedPtr mat = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols);
1515  (*mat) = df[2*dir][0]*stdiprod0 + df[2*dir+1][0]*stdiprod1;
1516 
1517  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
1518  }
1519  }
1520  break;
1521 
1522  case StdRegions::eInvHybridDGHelmholtz:
1523  {
1524  NekDouble one = 1.0;
1525 
1526  MatrixKey hkey(StdRegions::eHybridDGHelmholtz, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1527 
1528  DNekMatSharedPtr mat = GenMatrix(hkey);
1529 
1530  mat->Invert();
1531  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1532  }
1533  break;
1534  case StdRegions::ePreconLinearSpace:
1535  {
1536  NekDouble one = 1.0;
1537  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1538  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1539  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1540  DNekMatSharedPtr R=BuildVertexMatrix(A);
1541 
1542  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1543  }
1544  break;
1545  default:
1546  {
1547  NekDouble one = 1.0;
1548  DNekMatSharedPtr mat = GenMatrix(mkey);
1549 
1550  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1551  }
1552  break;
1553  }
1554 
1555  return returnval;
1556  }
1557 
1558 
1559  DNekScalBlkMatSharedPtr TriExp::CreateStaticCondMatrix(const MatrixKey &mkey)
1560  {
1561  DNekScalBlkMatSharedPtr returnval;
1562  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
1563 
1564  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
1565 
1566  // set up block matrix system
1567  unsigned int nbdry = NumBndryCoeffs();
1568  unsigned int nint = (unsigned int)(m_ncoeffs - nbdry);
1569  unsigned int exp_size[] = {nbdry,nint};
1570  unsigned int nblks = 2;
1571  returnval = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nblks,nblks,exp_size,exp_size); //Really need a constructor which takes Arrays
1572  NekDouble factor = 1.0;
1573 
1574  switch(mkey.GetMatrixType())
1575  {
1576  // this can only use stdregions statically condensed system for mass matrix
1577  case StdRegions::eMass:
1578  if((m_metricinfo->GetGtype() == SpatialDomains::eDeformed)||(mkey.GetNVarCoeff()))
1579  {
1580  factor = 1.0;
1581  goto UseLocRegionsMatrix;
1582  }
1583  else
1584  {
1585  factor = (m_metricinfo->GetJac(ptsKeys))[0];
1586  goto UseStdRegionsMatrix;
1587  }
1588  break;
1589  default: // use Deformed case for both regular and deformed geometries
1590  factor = 1.0;
1591  goto UseLocRegionsMatrix;
1592  break;
1593  UseStdRegionsMatrix:
1594  {
1595  NekDouble invfactor = 1.0/factor;
1596  NekDouble one = 1.0;
1597  DNekBlkMatSharedPtr mat = GetStdStaticCondMatrix(mkey);
1598  DNekScalMatSharedPtr Atmp;
1599  DNekMatSharedPtr Asubmat;
1600 
1601  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(0,0)));
1602  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,Asubmat = mat->GetBlock(0,1)));
1603  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(1,0)));
1604  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,Asubmat = mat->GetBlock(1,1)));
1605  }
1606  break;
1607 
1608  UseLocRegionsMatrix:
1609  {
1610  int i,j;
1611  NekDouble invfactor = 1.0/factor;
1612  NekDouble one = 1.0;
1613 
1614  DNekScalMat &mat = *GetLocMatrix(mkey);
1615 
1616  DNekMatSharedPtr A = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nbdry);
1617  DNekMatSharedPtr B = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nint);
1618  DNekMatSharedPtr C = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nbdry);
1619  DNekMatSharedPtr D = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nint);
1620 
1621  Array<OneD,unsigned int> bmap(nbdry);
1622  Array<OneD,unsigned int> imap(nint);
1623  GetBoundaryMap(bmap);
1624  GetInteriorMap(imap);
1625 
1626  for(i = 0; i < nbdry; ++i)
1627  {
1628  for(j = 0; j < nbdry; ++j)
1629  {
1630  (*A)(i,j) = mat(bmap[i],bmap[j]);
1631  }
1632 
1633  for(j = 0; j < nint; ++j)
1634  {
1635  (*B)(i,j) = mat(bmap[i],imap[j]);
1636  }
1637  }
1638 
1639  for(i = 0; i < nint; ++i)
1640  {
1641  for(j = 0; j < nbdry; ++j)
1642  {
1643  (*C)(i,j) = mat(imap[i],bmap[j]);
1644  }
1645 
1646  for(j = 0; j < nint; ++j)
1647  {
1648  (*D)(i,j) = mat(imap[i],imap[j]);
1649  }
1650  }
1651 
1652  // Calculate static condensed system
1653  if(nint)
1654  {
1655  D->Invert();
1656  (*B) = (*B)*(*D);
1657  (*A) = (*A) - (*B)*(*C);
1658  }
1659 
1660  DNekScalMatSharedPtr Atmp;
1661 
1662  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,A));
1663  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,B));
1664  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,C));
1665  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,D));
1666 
1667  }
1668  }
1669 
1670  return returnval;
1671  }
1672 
1673 
1674  DNekScalMatSharedPtr TriExp::v_GetLocMatrix(const MatrixKey &mkey)
1675  {
1676  return m_matrixManager[mkey];
1677  }
1678 
1679 
1680  DNekScalBlkMatSharedPtr TriExp::v_GetLocStaticCondMatrix(const MatrixKey &mkey)
1681  {
1682  return m_staticCondMatrixManager[mkey];
1683  }
1684 
1685  void TriExp::v_DropLocStaticCondMatrix(const MatrixKey &mkey)
1686  {
1687  m_staticCondMatrixManager.DeleteObject(mkey);
1688  }
1689 
1690 
1691 
1692  void TriExp::v_MassMatrixOp(const Array<OneD, const NekDouble> &inarray,
1693  Array<OneD,NekDouble> &outarray,
1694  const StdRegions::StdMatrixKey &mkey)
1695  {
1696  StdExpansion::MassMatrixOp_MatFree(inarray,outarray,mkey);
1697  }
1698 
1699 
1700  void TriExp::v_LaplacianMatrixOp(const Array<OneD, const NekDouble> &inarray,
1701  Array<OneD,NekDouble> &outarray,
1702  const StdRegions::StdMatrixKey &mkey)
1703  {
1704  TriExp::LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
1705  }
1706 
1707 
1708  void TriExp::v_LaplacianMatrixOp(const int k1, const int k2,
1709  const Array<OneD, const NekDouble> &inarray,
1710  Array<OneD,NekDouble> &outarray,
1711  const StdRegions::StdMatrixKey &mkey)
1712  {
1713  StdExpansion::LaplacianMatrixOp_MatFree(k1,k2,inarray,outarray,mkey);
1714  }
1715 
1716 
1717  void TriExp::v_WeakDerivMatrixOp(const int i,
1718  const Array<OneD, const NekDouble> &inarray,
1719  Array<OneD,NekDouble> &outarray,
1720  const StdRegions::StdMatrixKey &mkey)
1721  {
1722  StdExpansion::WeakDerivMatrixOp_MatFree(i,inarray,outarray,mkey);
1723  }
1724 
1725 
1726  void TriExp::v_WeakDirectionalDerivMatrixOp(const Array<OneD, const NekDouble> &inarray,
1727  Array<OneD,NekDouble> &outarray,
1728  const StdRegions::StdMatrixKey &mkey)
1729  {
1730  StdExpansion::WeakDirectionalDerivMatrixOp_MatFree(inarray,outarray,mkey);
1731  }
1732 
1733 
1734  void TriExp::v_MassLevelCurvatureMatrixOp(const Array<OneD, const NekDouble> &inarray,
1735  Array<OneD,NekDouble> &outarray,
1736  const StdRegions::StdMatrixKey &mkey)
1737  {
1738  StdExpansion::MassLevelCurvatureMatrixOp_MatFree(inarray,outarray,mkey);
1739  }
1740 
1741 
1742  void TriExp::v_HelmholtzMatrixOp(const Array<OneD, const NekDouble> &inarray,
1743  Array<OneD,NekDouble> &outarray,
1744  const StdRegions::StdMatrixKey &mkey)
1745  {
1746  TriExp::HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
1747  }
1748 
1749 
1750  void TriExp::v_GeneralMatrixOp_MatOp(const Array<OneD, const NekDouble> &inarray,
1751  Array<OneD,NekDouble> &outarray,
1752  const StdRegions::StdMatrixKey &mkey)
1753  {
1754  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1755 
1756  if(inarray.get() == outarray.get())
1757  {
1758  Array<OneD,NekDouble> tmp(m_ncoeffs);
1759  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
1760 
1761  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1762  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
1763  }
1764  else
1765  {
1766  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1767  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
1768  }
1769  }
1770 
1771 
1772  void TriExp::v_LaplacianMatrixOp_MatFree_Kernel(
1773  const Array<OneD, const NekDouble> &inarray,
1774  Array<OneD, NekDouble> &outarray,
1775  Array<OneD, NekDouble> &wsp)
1776  {
1777  if (m_metrics.count(eMetricLaplacian00) == 0)
1778  {
1779  ComputeLaplacianMetric();
1780  }
1781 
1782  int nquad0 = m_base[0]->GetNumPoints();
1783  int nquad1 = m_base[1]->GetNumPoints();
1784  int nqtot = nquad0*nquad1;
1785  int nmodes0 = m_base[0]->GetNumModes();
1786  int nmodes1 = m_base[1]->GetNumModes();
1787  int wspsize = max(max(max(nqtot,m_ncoeffs),nquad1*nmodes0),nquad0*nmodes1);
1788 
1789  ASSERTL1(wsp.num_elements() >= 3*wspsize,
1790  "Workspace is of insufficient size.");
1791 
1792  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
1793  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
1794  const Array<OneD, const NekDouble>& dbase0 = m_base[0]->GetDbdata();
1795  const Array<OneD, const NekDouble>& dbase1 = m_base[1]->GetDbdata();
1796  const Array<OneD, const NekDouble>& metric00 = m_metrics[eMetricLaplacian00];
1797  const Array<OneD, const NekDouble>& metric01 = m_metrics[eMetricLaplacian01];
1798  const Array<OneD, const NekDouble>& metric11 = m_metrics[eMetricLaplacian11];
1799 
1800  // Allocate temporary storage
1801  Array<OneD,NekDouble> wsp0(wsp);
1802  Array<OneD,NekDouble> wsp1(wsp+wspsize);
1803  Array<OneD,NekDouble> wsp2(wsp+2*wspsize);
1804 
1805  StdExpansion2D::PhysTensorDeriv(inarray,wsp1,wsp2);
1806 
1807  // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1808  // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1809  // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
1810  // especially for this purpose
1811  Vmath::Vvtvvtp(nqtot,&metric00[0],1,&wsp1[0],1,&metric01[0],1,&wsp2[0],1,&wsp0[0],1);
1812  Vmath::Vvtvvtp(nqtot,&metric01[0],1,&wsp1[0],1,&metric11[0],1,&wsp2[0],1,&wsp2[0],1);
1813 
1814  // outarray = m = (D_xi1 * B)^T * k
1815  // wsp1 = n = (D_xi2 * B)^T * l
1816  IProductWRTBase_SumFacKernel(dbase0,base1,wsp0,outarray,wsp1);
1817  IProductWRTBase_SumFacKernel(base0,dbase1,wsp2,wsp1, wsp0);
1818 
1819  // outarray = outarray + wsp1
1820  // = L * u_hat
1821  Vmath::Vadd(m_ncoeffs,wsp1.get(),1,outarray.get(),1,outarray.get(),1);
1822  }
1823 
1824 
1825  void TriExp::v_ComputeLaplacianMetric()
1826  {
1827  if (m_metrics.count(eMetricQuadrature) == 0)
1828  {
1829  ComputeQuadratureMetric();
1830  }
1831 
1832  unsigned int i, j;
1833  const SpatialDomains::GeomType type = m_metricinfo->GetGtype();
1834  const unsigned int nqtot = GetTotPoints();
1835  const unsigned int dim = 2;
1836  const MetricType m[3][3] = { {eMetricLaplacian00, eMetricLaplacian01, eMetricLaplacian02},
1837  {eMetricLaplacian01, eMetricLaplacian11, eMetricLaplacian12},
1838  {eMetricLaplacian02, eMetricLaplacian12, eMetricLaplacian22}
1839  };
1840 
1841  Array<OneD, NekDouble> dEta_dXi[2] = {Array<OneD, NekDouble>(nqtot,1.0),
1842  Array<OneD, NekDouble>(nqtot,1.0)};
1843 
1844  for (i = 0; i < dim; ++i)
1845  {
1846  for (j = i; j < dim; ++j)
1847  {
1848  m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
1849  }
1850  }
1851 
1852  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
1853  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
1854  const unsigned int nquad0 = m_base[0]->GetNumPoints();
1855  const unsigned int nquad1 = m_base[1]->GetNumPoints();
1856  const Array<TwoD, const NekDouble>& df =
1857  m_metricinfo->GetDerivFactors(GetPointsKeys());
1858 
1859  for(i = 0; i < nquad1; i++)
1860  {
1861  Blas::Dscal(nquad0,2.0/(1-z1[i]),&dEta_dXi[0][0]+i*nquad0,1);
1862  Blas::Dscal(nquad0,2.0/(1-z1[i]),&dEta_dXi[1][0]+i*nquad0,1);
1863  }
1864  for(i = 0; i < nquad0; i++)
1865  {
1866  Blas::Dscal(nquad1,0.5*(1+z0[i]),&dEta_dXi[1][0]+i,nquad0);
1867  }
1868 
1869  Array<OneD, NekDouble> tmp(nqtot);
1870  if((type == SpatialDomains::eRegular ||
1871  type == SpatialDomains::eMovingRegular))
1872  {
1873  Vmath::Smul (nqtot,df[0][0],&dEta_dXi[0][0],1,&tmp[0],1);
1874  Vmath::Svtvp(nqtot,df[1][0],&dEta_dXi[1][0],1,&tmp[0],1,&tmp[0],1);
1875 
1876  Vmath::Vmul (nqtot,&tmp[0],1, &tmp[0],1,&m_metrics[eMetricLaplacian00][0],1);
1877  Vmath::Smul (nqtot,df[1][0],&tmp[0],1,&m_metrics[eMetricLaplacian01][0],1);
1878 
1879 
1880  Vmath::Smul (nqtot,df[2][0],&dEta_dXi[0][0],1,&tmp[0],1);
1881  Vmath::Svtvp(nqtot,df[3][0],&dEta_dXi[1][0],1,&tmp[0],1,&tmp[0],1);
1882 
1883  Vmath::Vvtvp(nqtot,&tmp[0],1, &tmp[0],1,&m_metrics[eMetricLaplacian00][0],1,&m_metrics[eMetricLaplacian00][0],1);
1884  Vmath::Svtvp(nqtot,df[3][0],&tmp[0],1,&m_metrics[eMetricLaplacian01][0],1,&m_metrics[eMetricLaplacian01][0],1);
1885 
1886  if(GetCoordim() == 3)
1887  {
1888  Vmath::Smul (nqtot,df[4][0],&dEta_dXi[0][0],1,&tmp[0],1);
1889  Vmath::Svtvp(nqtot,df[5][0],&dEta_dXi[1][0],1,&tmp[0],1,&tmp[0],1);
1890 
1891  Vmath::Vvtvp(nqtot,&tmp[0],1, &tmp[0],1,&m_metrics[eMetricLaplacian00][0],1,&m_metrics[eMetricLaplacian00][0],1);
1892  Vmath::Svtvp(nqtot,df[5][0],&tmp[0],1,&m_metrics[eMetricLaplacian01][0],1,&m_metrics[eMetricLaplacian01][0],1);
1893  }
1894 
1895  NekDouble g2 = df[1][0]*df[1][0] + df[3][0]*df[3][0];
1896  if(GetCoordim() == 3)
1897  {
1898  g2 += df[5][0]*df[5][0];
1899  }
1900  Vmath::Fill(nqtot,g2,&m_metrics[eMetricLaplacian11][0],1);
1901  }
1902  else
1903  {
1904 
1905  Vmath::Vmul (nqtot,&df[0][0],1,&dEta_dXi[0][0],1,&tmp[0],1);
1906  Vmath::Vvtvp(nqtot,&df[1][0],1,&dEta_dXi[1][0],1,&tmp[0],1,&tmp[0],1);
1907 
1908  Vmath::Vmul (nqtot,&tmp[0], 1,&tmp[0], 1,&m_metrics[eMetricLaplacian00][0],1);
1909  Vmath::Vmul (nqtot,&df[1][0],1,&tmp[0], 1,&m_metrics[eMetricLaplacian01][0],1);
1910  Vmath::Vmul (nqtot,&df[1][0],1,&df[1][0],1,&m_metrics[eMetricLaplacian11][0],1);
1911 
1912 
1913  Vmath::Vmul (nqtot,&df[2][0],1,&dEta_dXi[0][0],1,&tmp[0],1);
1914  Vmath::Vvtvp(nqtot,&df[3][0],1,&dEta_dXi[1][0],1,&tmp[0],1,&tmp[0],1);
1915 
1916  Vmath::Vvtvp(nqtot,&tmp[0], 1,&tmp[0], 1,&m_metrics[eMetricLaplacian00][0],1,&m_metrics[eMetricLaplacian00][0],1);
1917  Vmath::Vvtvp(nqtot,&df[3][0],1,&tmp[0], 1,&m_metrics[eMetricLaplacian01][0],1,&m_metrics[eMetricLaplacian01][0],1);
1918  Vmath::Vvtvp(nqtot,&df[3][0],1,&df[3][0],1,&m_metrics[eMetricLaplacian11][0],1,&m_metrics[eMetricLaplacian11][0],1);
1919 
1920  if(GetCoordim() == 3)
1921  {
1922  Vmath::Vmul (nqtot,&df[4][0],1,&dEta_dXi[0][0],1,&tmp[0],1);
1923  Vmath::Vvtvp(nqtot,&df[5][0],1,&dEta_dXi[1][0],1,&tmp[0],1,&tmp[0],1);
1924 
1925  Vmath::Vvtvp(nqtot,&tmp[0], 1,&tmp[0], 1,&m_metrics[eMetricLaplacian00][0],1,&m_metrics[eMetricLaplacian00][0],1);
1926  Vmath::Vvtvp(nqtot,&df[5][0],1,&tmp[0], 1,&m_metrics[eMetricLaplacian01][0],1,&m_metrics[eMetricLaplacian01][0],1);
1927  Vmath::Vvtvp(nqtot,&df[5][0],1,&df[5][0],1,&m_metrics[eMetricLaplacian11][0],1,&m_metrics[eMetricLaplacian11][0],1);
1928  }
1929  }
1930 
1931  for (unsigned int i = 0; i < dim; ++i)
1932  {
1933  for (unsigned int j = i; j < dim; ++j)
1934  {
1935  MultiplyByQuadratureMetric(m_metrics[m[i][j]],
1936  m_metrics[m[i][j]]);
1937 
1938  }
1939  }
1940  }
1941 
1942  /**
1943  * Function is used to compute exactly the advective numerical flux on
1944  * theinterface of two elements with different expansions, hence an
1945  * appropriate number of Gauss points has to be used. The number of
1946  * Gauss points has to be equal to the number used by the highest
1947  * polynomial degree of the two adjacent elements. Furthermore, this
1948  * function is used to compute the sensor value in each element.
1949  *
1950  * @param numMin Is the reduced polynomial order
1951  * @param inarray Input array of coefficients
1952  * @param dumpVar Output array of reduced coefficients.
1953  */
1954  void TriExp::v_ReduceOrderCoeffs(
1955  int numMin,
1956  const Array<OneD, const NekDouble> &inarray,
1957  Array<OneD, NekDouble> &outarray)
1958  {
1959  int n_coeffs = inarray.num_elements();
1960  int nquad0 = m_base[0]->GetNumPoints();
1961  int nquad1 = m_base[1]->GetNumPoints();
1962  int nqtot = nquad0*nquad1;
1963  int nmodes0 = m_base[0]->GetNumModes();
1964  int nmodes1 = m_base[1]->GetNumModes();
1965  int numMin2 = nmodes0, i;
1966 
1967  Array<OneD, NekDouble> coeff(n_coeffs,0.0);
1968  Array<OneD, NekDouble> phys_tmp(nqtot,0.0);
1969  Array<OneD, NekDouble> tmp, tmp2;
1970 
1971  const LibUtilities::PointsKey Pkey0 = m_base[0]->GetPointsKey();
1972  const LibUtilities::PointsKey Pkey1 = m_base[1]->GetPointsKey();
1973 
1974  LibUtilities::BasisKey b0(
1975  m_base[0]->GetBasisType(), nmodes0, Pkey0);
1976  LibUtilities::BasisKey b1(
1977  m_base[1]->GetBasisType(), nmodes1, Pkey1);
1978  LibUtilities::BasisKey bortho0(
1979  LibUtilities::eOrtho_A, nmodes0, Pkey0);
1980  LibUtilities::BasisKey bortho1(
1981  LibUtilities::eOrtho_B, nmodes1, Pkey1);
1982 
1983  // Check if it is also possible to use the same InterCoeff routine
1984  // which is also used for Quadrilateral and Hexagonal shaped
1985  // elements
1986 
1987  // For now, set up the used basis on the standard element to
1988  // calculate the phys values, set up the orthogonal basis to do a
1989  // forward transform, to obtain the coefficients in orthogonal
1990  // coefficient space
1991  StdRegions::StdTriExpSharedPtr m_OrthoTriExp;
1992  StdRegions::StdTriExpSharedPtr m_TriExp;
1993 
1994  m_TriExp = MemoryManager<StdRegions::StdTriExp>
1995  ::AllocateSharedPtr(b0, b1);
1996  m_OrthoTriExp = MemoryManager<StdRegions::StdTriExp>
1997  ::AllocateSharedPtr(bortho0, bortho1);
1998 
1999  m_TriExp ->BwdTrans(inarray,phys_tmp);
2000  m_OrthoTriExp->FwdTrans(phys_tmp, coeff);
2001 
2002  for (i = 0; i < n_coeffs; i++)
2003  {
2004  if (i == numMin)
2005  {
2006  coeff[i] = 0.0;
2007  numMin += numMin2 - 1;
2008  numMin2 -= 1.0;
2009  }
2010  }
2011 
2012  m_OrthoTriExp->BwdTrans(coeff,phys_tmp);
2013  m_TriExp ->FwdTrans(phys_tmp, outarray);
2014  }
2015 
2016  void TriExp::v_SVVLaplacianFilter(
2017  Array<OneD, NekDouble> &array,
2018  const StdRegions::StdMatrixKey &mkey)
2019  {
2020  int nq = GetTotPoints();
2021 
2022  // Calculate sqrt of the Jacobian
2023  Array<OneD, const NekDouble> jac =
2024  m_metricinfo->GetJac(GetPointsKeys());
2025  Array<OneD, NekDouble> sqrt_jac(nq);
2026  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
2027  {
2028  Vmath::Vsqrt(nq,jac,1,sqrt_jac,1);
2029  }
2030  else
2031  {
2032  Vmath::Fill(nq,sqrt(jac[0]),sqrt_jac,1);
2033  }
2034 
2035  // Multiply array by sqrt(Jac)
2036  Vmath::Vmul(nq,sqrt_jac,1,array,1,array,1);
2037 
2038  // Apply std region filter
2039  StdTriExp::v_SVVLaplacianFilter( array, mkey);
2040 
2041  // Divide by sqrt(Jac)
2042  Vmath::Vdiv(nq,array,1,sqrt_jac,1,array,1);
2043  }
2044 
2045  }
2046 }
2047 
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:250
ASSERTL2
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:274
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:16
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:60
Nektar::LibUtilities::PointsKey::GetNumPoints
unsigned int GetNumPoints() const
Definition: Points.h:107
Nektar::LocalRegions::Expansion2D::GetLeftAdjacentElementExp
Expansion3DSharedPtr GetLeftAdjacentElementExp() const
Definition: Expansion2D.h:245
Nektar::LocalRegions::Expansion2D::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: Expansion2D.cpp:713
Nektar::LocalRegions::Expansion2D::GetGeom2D
SpatialDomains::Geometry2DSharedPtr GetGeom2D() const
Definition: Expansion2D.h:291
Nektar::LocalRegions::Expansion2D::m_edgeNormals
std::map< int, StdRegions::NormalVector > m_edgeNormals
Definition: Expansion2D.h:134
Nektar::LocalRegions::Expansion::GetGeom
SpatialDomains::GeometrySharedPtr GetGeom() const
Definition: Expansion.cpp:167
Nektar::LocalRegions::Expansion::m_geom
SpatialDomains::GeometrySharedPtr m_geom
Definition: Expansion.h:127
Nektar::LocalRegions::Expansion::BuildVertexMatrix
DNekMatSharedPtr BuildVertexMatrix(const DNekScalMatSharedPtr &r_bnd)
Definition: Expansion.cpp:98
Nektar::LocalRegions::Expansion::ComputeGmatcdotMF
void ComputeGmatcdotMF(const Array< TwoD, const NekDouble > &df, const Array< OneD, const NekDouble > &direction, Array< OneD, Array< OneD, NekDouble > > &dfdir)
Definition: Expansion.cpp:312
Nektar::LocalRegions::Expansion::m_metricinfo
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
Definition: Expansion.h:128
Nektar::LocalRegions::Expansion::v_GetCoords
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
Definition: Expansion.cpp:231
Nektar::LocalRegions::Expansion::GetLocMatrix
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition: Expansion.cpp:85
Nektar::LocalRegions::TriExp::m_matrixManager
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
Definition: TriExp.h:296
Nektar::LocalRegions::TriExp::v_PhysDirectionalDeriv
virtual void v_PhysDirectionalDeriv(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &direction, Array< OneD, NekDouble > &out)
Physical derivative along a direction vector.
Definition: TriExp.cpp:196
Nektar::LocalRegions::TriExp::v_NormVectorIProductWRTBase
virtual void v_NormVectorIProductWRTBase(const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, const Array< OneD, const NekDouble > &Fz, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:655
Nektar::LocalRegions::TriExp::v_GetCoord
virtual void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords)
Definition: TriExp.cpp:710
Nektar::LocalRegions::TriExp::v_GetEdgePhysVals
virtual void v_GetEdgePhysVals(const int edge, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Extract the physical values along edge edge from inarray into outarray following the local edge orien...
Definition: TriExp.cpp:771
Nektar::LocalRegions::TriExp::v_LaplacianMatrixOp
virtual void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1700
Nektar::LocalRegions::TriExp::v_ExtractDataToCoeffs
virtual void v_ExtractDataToCoeffs(const NekDouble *data, const std::vector< unsigned int > &nummodes, const int mode_offset, NekDouble *coeffs, std::vector< LibUtilities::BasisType > &fromType)
Definition: TriExp.cpp:1111
Nektar::LocalRegions::TriExp::v_WeakDerivMatrixOp
virtual void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1717
Nektar::LocalRegions::TriExp::v_MassMatrixOp
virtual void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1692
Nektar::LocalRegions::TriExp::v_GetLocMatrix
virtual DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey)
Definition: TriExp.cpp:1674
Nektar::LocalRegions::TriExp::v_ReduceOrderCoeffs
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:1954
Nektar::LocalRegions::TriExp::v_GetBasis
virtual const LibUtilities::BasisSharedPtr & v_GetBasis(int dir) const
Definition: TriExp.cpp:1157
Nektar::LocalRegions::TriExp::v_IProductWRTDerivBase
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:414
Nektar::LocalRegions::TriExp::v_GetLocStaticCondMatrix
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey)
Definition: TriExp.cpp:1680
Nektar::LocalRegions::TriExp::v_PhysDeriv
virtual void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
Calculate the derivative of the physical points.
Definition: TriExp.cpp:108
Nektar::LocalRegions::TriExp::m_staticCondMatrixManager
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
Definition: TriExp.h:297
Nektar::LocalRegions::TriExp::v_GetCoordim
virtual int v_GetCoordim()
Definition: TriExp.cpp:1105
Nektar::LocalRegions::TriExp::v_GetNumPoints
virtual int v_GetNumPoints(const int dir) const
Definition: TriExp.cpp:1164
Nektar::LocalRegions::TriExp::v_IProductWRTBase_SumFac
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Definition: TriExp.cpp:422
Nektar::LocalRegions::TriExp::v_GetEdgeInterpVals
virtual void v_GetEdgeInterpVals(const int edge, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:880
Nektar::LocalRegions::TriExp::v_GetStdExp
virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const
Definition: TriExp.cpp:691
Nektar::LocalRegions::TriExp::v_DropLocStaticCondMatrix
void v_DropLocStaticCondMatrix(const MatrixKey &mkey)
Definition: TriExp.cpp:1685
Nektar::LocalRegions::TriExp::v_Integral
virtual NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
Integrates the specified function over the domain.
Definition: TriExp.cpp:85
Nektar::LocalRegions::TriExp::v_MassLevelCurvatureMatrixOp
virtual void v_MassLevelCurvatureMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1734
Nektar::LocalRegions::TriExp::v_ComputeLaplacianMetric
virtual void v_ComputeLaplacianMetric()
Definition: TriExp.cpp:1825
Nektar::LocalRegions::TriExp::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1170
Nektar::LocalRegions::TriExp::v_SVVLaplacianFilter
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:2016
Nektar::LocalRegions::TriExp::v_GetEorient
virtual StdRegions::Orientation v_GetEorient(int edge)
Definition: TriExp.cpp:1152
Nektar::LocalRegions::TriExp::v_GetEdgeQFactors
virtual void v_GetEdgeQFactors(const int edge, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:889
Nektar::LocalRegions::TriExp::v_LaplacianMatrixOp_MatFree_Kernel
virtual void v_LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp)
Definition: TriExp.cpp:1772
Nektar::LocalRegions::TriExp::v_CreateStdMatrix
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1193
Nektar::LocalRegions::TriExp::v_IProductWRTBase_MatOp
virtual void v_IProductWRTBase_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:448
Nektar::LocalRegions::TriExp::v_FwdTrans
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Transform a given function from physical quadrature space to coefficient space.
Definition: TriExp.cpp:267
Nektar::LocalRegions::TriExp::v_IProductWRTBase
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Calculate the inner product of inarray with respect to the basis B=base0[p]*base1[pq] and put into ou...
Definition: TriExp.cpp:407
Nektar::LocalRegions::TriExp::v_PhysEvaluate
virtual NekDouble v_PhysEvaluate(const Array< OneD, const NekDouble > &coord, const Array< OneD, const NekDouble > &physvals)
This function evaluates the expansion at a single (arbitrary) point of the domain.
Definition: TriExp.cpp:749
Nektar::LocalRegions::TriExp::v_StdPhysEvaluate
virtual NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals)
Definition: TriExp.cpp:741
Nektar::LocalRegions::TriExp::v_ComputeEdgeNormal
virtual void v_ComputeEdgeNormal(const int edge)
Definition: TriExp.cpp:939
Nektar::LocalRegions::TriExp::CreateMatrix
virtual DNekScalMatSharedPtr CreateMatrix(const MatrixKey &mkey)
Definition: TriExp.cpp:1204
Nektar::LocalRegions::TriExp::v_GetCoords
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
Definition: TriExp.cpp:727
Nektar::LocalRegions::TriExp::v_HelmholtzMatrixOp
virtual void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1742
Nektar::LocalRegions::TriExp::v_GetTracePhysVals
virtual void v_GetTracePhysVals(const int edge, const StdRegions::StdExpansionSharedPtr &EdgeExp, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, StdRegions::Orientation orient)
Definition: TriExp.cpp:760
Nektar::LocalRegions::TriExp::v_IProductWRTDirectionalDerivBase
virtual void v_IProductWRTDirectionalDerivBase(const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:571
Nektar::LocalRegions::TriExp::v_FwdTrans_BndConstrained
virtual void v_FwdTrans_BndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:285
Nektar::LocalRegions::TriExp::v_GetLinStdExp
virtual StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const
Definition: TriExp.cpp:699
Nektar::LocalRegions::TriExp::v_IProductWRTDirectionalDerivBase_SumFac
virtual void v_IProductWRTDirectionalDerivBase_SumFac(const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Directinoal Derivative in the modal space in the dir direction of varcoeffs.
Definition: TriExp.cpp:585
Nektar::LocalRegions::TriExp::v_GetEdgePhysMap
virtual void v_GetEdgePhysMap(const int edge, Array< OneD, int > &outarray)
Definition: TriExp.cpp:900
Nektar::LocalRegions::TriExp::CreateStaticCondMatrix
virtual DNekScalBlkMatSharedPtr CreateStaticCondMatrix(const MatrixKey &mkey)
Definition: TriExp.cpp:1559
Nektar::LocalRegions::TriExp::v_GeneralMatrixOp_MatOp
virtual void v_GeneralMatrixOp_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1750
Nektar::LocalRegions::TriExp::v_WeakDirectionalDerivMatrixOp
virtual void v_WeakDirectionalDerivMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1726
Nektar::LocalRegions::TriExp::v_IProductWRTDerivBase_MatOp
virtual void v_IProductWRTDerivBase_MatOp(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:531
Nektar::LocalRegions::TriExp::v_IProductWRTDerivBase_SumFac
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:461
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:84
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:161
Nektar::StdRegions::StdExpansion2D::IProductWRTBase_SumFacKernel
void IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0=true, bool doCheckCollDir1=true)
Definition: StdExpansion2D.cpp:199
Nektar::StdRegions::StdExpansion::GetBoundaryMap
void GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:812
Nektar::StdRegions::StdExpansion::GetTotPoints
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:140
Nektar::StdRegions::StdExpansion::GetBasisType
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
Definition: StdExpansion.h:164
Nektar::StdRegions::StdExpansion::HelmholtzMatrixOp_MatFree
void HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:1550
Nektar::StdRegions::StdExpansion::GetStdMatrix
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:714
Nektar::StdRegions::StdExpansion::GetPointsKeys
const LibUtilities::PointsKeyVector GetPointsKeys() const
Definition: StdExpansion.h:1330
Nektar::StdRegions::StdExpansion::GetStdStaticCondMatrix
DNekBlkMatSharedPtr GetStdStaticCondMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:719
Nektar::StdRegions::StdExpansion::MassMatrixOp
void MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:974
Nektar::StdRegions::StdExpansion::GetLocStaticCondMatrix
DNekScalBlkMatSharedPtr GetLocStaticCondMatrix(const LocalRegions::MatrixKey &mkey)
Definition: StdExpansion.h:761
Nektar::StdRegions::StdExpansion::IProductWRTDerivBase_SumFac
void IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:1474
Nektar::StdRegions::StdExpansion::GetEorient
StdRegions::Orientation GetEorient(int edge)
Definition: StdExpansion.h:776
Nektar::StdRegions::StdExpansion::NormVectorIProductWRTBase
void NormVectorIProductWRTBase(const Array< OneD, const NekDouble > &Fx, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:741
Nektar::StdRegions::StdExpansion::IProductWRTBase
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
this function calculates the inner product of a given function f with the different modes of the expa...
Definition: StdExpansion.h:634
Nektar::StdRegions::StdExpansion::LaplacianMatrixOp_MatFree
void LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:1508
Nektar::StdRegions::StdExpansion::DetShapeType
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:469
Nektar::StdRegions::StdExpansion::GetEdgeToElementMap
void GetEdgeToElementMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, int P=-1)
Definition: StdExpansion.h:849
Nektar::StdRegions::StdExpansion::GetInteriorMap
void GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:817
Nektar::StdRegions::StdExpansion::GenMatrix
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:1060
Nektar::StdRegions::StdExpansion::GetPointsType
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:215
Nektar::StdRegions::StdExpansion::GetNumPoints
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
Definition: StdExpansion.h:228
Nektar::StdRegions::StdExpansion::IProductWRTDirectionalDerivBase_SumFac
void IProductWRTDirectionalDerivBase_SumFac(const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:1483
Nektar::StdRegions::StdExpansion::IProductWRTBase_SumFac
void IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Definition: StdExpansion.h:1421
Nektar::StdRegions::StdExpansion::MultiplyByQuadratureMetric
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:945
Nektar::StdRegions::StdExpansion::PhysDeriv
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
Definition: StdExpansion.h:1065
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition: StdExpansion.h:1429
Nektar::StdRegions::StdMatrixKey::GetConstFactors
const ConstFactorMap & GetConstFactors() const
Definition: StdMatrixKey.h:135
Nektar::StdRegions::StdMatrixKey::GetShapeType
LibUtilities::ShapeType GetShapeType() const
Definition: StdMatrixKey.h:86
Nektar::StdRegions::StdMatrixKey::GetVarCoeffs
const VarCoeffMap & GetVarCoeffs() const
Definition: StdMatrixKey.h:161
Nektar::StdRegions::StdMatrixKey::GetMatrixType
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:81
Nektar::StdRegions::StdMatrixKey::GetVarCoeff
const Array< OneD, const NekDouble > & GetVarCoeff(const StdRegions::VarCoeffType &coeff) const
Definition: StdMatrixKey.h:145
Nektar::StdRegions::StdMatrixKey::GetConstFactor
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:121
Nektar::StdRegions::StdMatrixKey::ConstFactorExists
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:130
Blas::Dgemv
static void Dgemv(const char &trans, const int &m, const int &n, const double &alpha, const double *a, const int &lda, const double *x, const int &incx, const double &beta, double *y, const int &incy)
BLAS level 2: Matrix vector multiply y = A x where A[m x n].
Definition: Blas.hpp:168
Blas::Dscal
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:125
Blas::Daxpy
static void Daxpy(const int &n, const double &alpha, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: y = alpha x plus y.
Definition: Blas.hpp:110
Nektar::LibUtilities::BasisSharedPtr
std::shared_ptr< Basis > BasisSharedPtr
Definition: FoundationsFwd.hpp:74
Nektar::LibUtilities::Interp1D
void Interp1D(const BasisKey &fbasis0, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, Array< OneD, NekDouble > &to)
this function interpolates a 1D function evaluated at the quadrature points of the basis fbasis0 to ...
Definition: Interp.cpp:53
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:246
Nektar::LibUtilities::eGaussLobattoLegendre
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:51
Nektar::LibUtilities::eModified_B
@ eModified_B
Principle Modified Functions .
Definition: BasisType.h:49
Nektar::LibUtilities::eOrtho_A
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:45
Nektar::LibUtilities::eOrtho_B
@ eOrtho_B
Principle Orthogonal Functions .
Definition: BasisType.h:46
Nektar::LibUtilities::eModified_A
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
Nektar::LocalRegions::SegExpSharedPtr
std::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:266
Nektar::SpatialDomains::GeomFactorsSharedPtr
std::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
Nektar::SpatialDomains::GeomType
GeomType
Indicates the type of element geometry.
Definition: SpatialDomains.hpp:52
Nektar::SpatialDomains::eRegular
@ eRegular
Geometry is straight-sided with constant geometric factors.
Definition: SpatialDomains.hpp:54
Nektar::SpatialDomains::eNoGeomType
@ eNoGeomType
No type defined.
Definition: SpatialDomains.hpp:53
Nektar::SpatialDomains::eMovingRegular
@ eMovingRegular
Currently unused.
Definition: SpatialDomains.hpp:57
Nektar::SpatialDomains::eDeformed
@ eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:56
Nektar::SpatialDomains::TriGeomSharedPtr
std::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
Nektar::StdRegions::StdExpansionSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1926
Nektar::StdRegions::VarCoeffMap
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
Definition: StdRegions.hpp:264
Nektar::StdRegions::StdTriExpSharedPtr
std::shared_ptr< StdTriExp > StdTriExpSharedPtr
Definition: StdTriExp.h:266
Nektar::StdRegions::NullConstFactorMap
static ConstFactorMap NullConstFactorMap
Definition: StdRegions.hpp:295
Nektar::rhs
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs
Definition: MatrixOperations.cpp:97
Nektar::DNekScalMatSharedPtr
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Definition: NekMatrixFwd.hpp:69
Nektar::DNekBlkMatSharedPtr
std::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
Definition: NekTypeDefs.hpp:71
Nektar::Transpose
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Definition: BlockMatrix.hpp:217
Nektar::DNekScalBlkMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:73
Nektar::NullNekDouble1DArray
static Array< OneD, NekDouble > NullNekDouble1DArray
Definition: SharedArray.hpp:787
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:69
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:411
Vmath::Svtsvtp
void Svtsvtp(int n, const T alpha, const T *x, int incx, const T beta, const T *y, int incy, T *z, int incz)
vvtvvtp (scalar times vector plus scalar times vector):
Definition: Vmath.cpp:594
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:488
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:445
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:216
Vmath::Sdiv
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
Definition: Vmath.cpp:274
Vmath::Vdiv
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:244
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:376
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1088
Vmath::Vvtvvtp
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
Definition: Vmath.cpp:540
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:346
Generated on Mon Aug 23 2021 06:57:51 for Nektar++ by doxygen 1.9.1