Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
QuadExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File QuadExp.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: Expansion for quadrilateral elements.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 #include <LibUtilities/Foundations/InterpCoeff.h>
36 #include <LocalRegions/QuadExp.h>
37 #include <LocalRegions/Expansion3D.h>
38 #include <LibUtilities/BasicUtils/VmathArray.hpp>
39 #include <LibUtilities/BasicUtils/Vmath.hpp>
40 #include <LibUtilities/Foundations/Interp.h>
41 #include <LocalRegions/SegExp.h>
42 
43 using namespace std;
44 
45 namespace Nektar
46 {
47  namespace LocalRegions
48  {
49  QuadExp::QuadExp(const LibUtilities::BasisKey &Ba,
50  const LibUtilities::BasisKey &Bb,
51  const SpatialDomains::QuadGeomSharedPtr &geom):
52  StdExpansion (Ba.GetNumModes()*Bb.GetNumModes(),2,Ba,Bb),
53  StdExpansion2D(Ba.GetNumModes()*Bb.GetNumModes(),Ba,Bb),
54  StdQuadExp (Ba,Bb),
55  Expansion (geom),
56  Expansion2D (geom),
57  m_matrixManager(
58  boost::bind(&QuadExp::CreateMatrix, this, _1),
59  std::string("QuadExpMatrix")),
60  m_staticCondMatrixManager(
61  boost::bind(&QuadExp::CreateStaticCondMatrix, this, _1),
62  std::string("QuadExpStaticCondMatrix"))
63  {
64  }
65 
66 
67  QuadExp::QuadExp(const QuadExp &T):
68  StdExpansion(T),
69  StdExpansion2D(T),
70  StdQuadExp(T),
71  Expansion (T),
72  Expansion2D (T),
73  m_matrixManager(T.m_matrixManager),
74  m_staticCondMatrixManager(T.m_staticCondMatrixManager)
75  {
76  }
77 
78 
79  QuadExp::~QuadExp()
80  {
81  }
82 
83 
84  NekDouble QuadExp::v_Integral(
85  const Array<OneD,
86  const NekDouble> &inarray)
87  {
88  int nquad0 = m_base[0]->GetNumPoints();
89  int nquad1 = m_base[1]->GetNumPoints();
90  Array<OneD, const NekDouble> jac = m_metricinfo->GetJac(GetPointsKeys());
91  NekDouble ival;
92  Array<OneD,NekDouble> tmp(nquad0*nquad1);
93 
94  // multiply inarray with Jacobian
95  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
96  {
97  Vmath::Vmul(nquad0*nquad1, jac, 1, inarray, 1, tmp, 1);
98  }
99  else
100  {
101  Vmath::Smul(nquad0*nquad1, jac[0], inarray, 1, tmp, 1);
102  }
103 
104  // call StdQuadExp version;
105  ival = StdQuadExp::v_Integral(tmp);
106  return ival;
107  }
108 
109 
110  void QuadExp::v_PhysDeriv(
111  const Array<OneD, const NekDouble> & inarray,
112  Array<OneD,NekDouble> &out_d0,
113  Array<OneD,NekDouble> &out_d1,
114  Array<OneD,NekDouble> &out_d2)
115  {
116  int nquad0 = m_base[0]->GetNumPoints();
117  int nquad1 = m_base[1]->GetNumPoints();
118  int nqtot = nquad0*nquad1;
119  const Array<TwoD, const NekDouble>& df = m_metricinfo->GetDerivFactors(GetPointsKeys());
120  Array<OneD,NekDouble> diff0(2*nqtot);
121  Array<OneD,NekDouble> diff1(diff0+nqtot);
122 
123  StdQuadExp::v_PhysDeriv(inarray, diff0, diff1);
124 
125  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
126  {
127  if (out_d0.num_elements())
128  {
129  Vmath::Vmul (nqtot, df[0], 1, diff0, 1, out_d0, 1);
130  Vmath::Vvtvp (nqtot, df[1], 1, diff1, 1, out_d0, 1,
131  out_d0,1);
132  }
133 
134  if(out_d1.num_elements())
135  {
136  Vmath::Vmul (nqtot,df[2],1,diff0,1, out_d1, 1);
137  Vmath::Vvtvp (nqtot,df[3],1,diff1,1, out_d1, 1, out_d1,1);
138  }
139 
140  if (out_d2.num_elements())
141  {
142  Vmath::Vmul (nqtot,df[4],1,diff0,1, out_d2, 1);
143  Vmath::Vvtvp (nqtot,df[5],1,diff1,1, out_d2, 1, out_d2,1);
144  }
145  }
146  else // regular geometry
147  {
148  if (out_d0.num_elements())
149  {
150  Vmath::Smul (nqtot, df[0][0], diff0, 1, out_d0, 1);
151  Blas::Daxpy (nqtot, df[1][0], diff1, 1, out_d0, 1);
152  }
153 
154  if (out_d1.num_elements())
155  {
156  Vmath::Smul (nqtot, df[2][0], diff0, 1, out_d1, 1);
157  Blas::Daxpy (nqtot, df[3][0], diff1, 1, out_d1, 1);
158  }
159 
160  if (out_d2.num_elements())
161  {
162  Vmath::Smul (nqtot, df[4][0], diff0, 1, out_d2, 1);
163  Blas::Daxpy (nqtot, df[5][0], diff1, 1, out_d2, 1);
164  }
165  }
166  }
167 
168 
169  void QuadExp::v_PhysDeriv(
170  const int dir,
171  const Array<OneD, const NekDouble>& inarray,
172  Array<OneD, NekDouble> &outarray)
173  {
174  switch (dir)
175  {
176  case 0:
177  {
178  v_PhysDeriv(inarray, outarray, NullNekDouble1DArray,
179  NullNekDouble1DArray);
180  }
181  break;
182  case 1:
183  {
184  v_PhysDeriv(inarray, NullNekDouble1DArray, outarray,
185  NullNekDouble1DArray);
186  }
187  break;
188  case 2:
189  {
190  v_PhysDeriv(inarray, NullNekDouble1DArray,
191  NullNekDouble1DArray, outarray);
192  }
193  break;
194  default:
195  {
196  ASSERTL1(false,"input dir is out of range");
197  }
198  break;
199  }
200  }
201 
202 
203  void QuadExp::v_PhysDirectionalDeriv(
204  const Array<OneD, const NekDouble> & inarray,
205  const Array<OneD, const NekDouble>& direction,
206  Array<OneD,NekDouble> &out)
207  {
208  int nquad0 = m_base[0]->GetNumPoints();
209  int nquad1 = m_base[1]->GetNumPoints();
210  int nqtot = nquad0*nquad1;
211 
212  const Array<TwoD, const NekDouble>& df = m_metricinfo->GetDerivFactors(GetPointsKeys());
213 
214  Array<OneD,NekDouble> diff0(2*nqtot);
215  Array<OneD,NekDouble> diff1(diff0+nqtot);
216 
217  StdQuadExp::v_PhysDeriv(inarray, diff0, diff1);
218 
219  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
220  {
221  Array<OneD, Array<OneD, NekDouble> > tangmat(2);
222 
223  // d/dx_v^s = v_x*ds/dx + v_y*ds/dy + v_z*dx/dz
224  for (int i=0; i< 2; ++i)
225  {
226  tangmat[i] = Array<OneD, NekDouble>(nqtot,0.0);
227  for (int k=0; k<(m_geom->GetCoordim()); ++k)
228  {
229  Vmath::Vvtvp(nqtot,
230  &df[2*k+i][0], 1,
231  &direction[k*nqtot], 1,
232  &tangmat[i][0], 1,
233  &tangmat[i][0], 1);
234  }
235  }
236 
237  /// D_v = d/dx_v^s + d/dx_v^r
238  if (out.num_elements())
239  {
240  Vmath::Vmul (nqtot,
241  &tangmat[0][0], 1,
242  &diff0[0], 1,
243  &out[0], 1);
244  Vmath::Vvtvp (nqtot,
245  &tangmat[1][0], 1,
246  &diff1[0], 1,
247  &out[0], 1,
248  &out[0], 1);
249  }
250 
251  }
252  else
253  {
254  ASSERTL1(m_metricinfo->GetGtype() ==
255  SpatialDomains::eDeformed,"Wrong route");
256  }
257  }
258 
259 
260  void QuadExp::v_FwdTrans(
261  const Array<OneD, const NekDouble> & inarray,
262  Array<OneD,NekDouble> &outarray)
263  {
264  if ((m_base[0]->Collocation())&&(m_base[1]->Collocation()))
265  {
266  Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
267  }
268  else
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 
286  void QuadExp::v_FwdTrans_BndConstrained(
287  const Array<OneD, const NekDouble>& inarray,
288  Array<OneD, NekDouble> &outarray)
289  {
290  if ((m_base[0]->Collocation())&&(m_base[1]->Collocation()))
291  {
292  Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
293  }
294  else
295  {
296  int i,j;
297  int npoints[2] = {m_base[0]->GetNumPoints(),
298  m_base[1]->GetNumPoints()};
299  int nmodes[2] = {m_base[0]->GetNumModes(),
300  m_base[1]->GetNumModes()};
301 
302  fill(outarray.get(), outarray.get()+m_ncoeffs, 0.0 );
303 
304  Array<OneD, NekDouble> physEdge[4];
305  Array<OneD, NekDouble> coeffEdge[4];
306  StdRegions::Orientation orient[4];
307  for (i = 0; i < 4; i++)
308  {
309  physEdge[i] = Array<OneD, NekDouble>(npoints[i%2]);
310  coeffEdge[i] = Array<OneD, NekDouble>(nmodes[i%2]);
311  orient[i] = GetEorient(i);
312  }
313 
314  for (i = 0; i < npoints[0]; i++)
315  {
316  physEdge[0][i] = inarray[i];
317  physEdge[2][i] = inarray[npoints[0]*npoints[1]-1-i];
318  }
319 
320  for (i = 0; i < npoints[1]; i++)
321  {
322  physEdge[1][i] =
323  inarray[npoints[0]-1+i*npoints[0]];
324  physEdge[3][i] =
325  inarray[(npoints[1]-1)*npoints[0]-i*npoints[0]];
326  }
327 
328  for (i = 0; i < 4; i++)
329  {
330  if ( orient[i] == StdRegions::eBackwards )
331  {
332  reverse((physEdge[i]).get(),
333  (physEdge[i]).get() + npoints[i%2] );
334  }
335  }
336 
337  SegExpSharedPtr segexp[4];
338  for (i = 0; i < 4; i++)
339  {
340  segexp[i] = MemoryManager<LocalRegions::SegExp>::
341  AllocateSharedPtr(
342  m_base[i%2]->GetBasisKey(),GetGeom2D()->GetEdge(i));
343  }
344 
345  Array<OneD, unsigned int> mapArray;
346  Array<OneD, int> signArray;
347  NekDouble sign;
348 
349  for (i = 0; i < 4; i++)
350  {
351  segexp[i%2]->FwdTrans_BndConstrained(
352  physEdge[i],coeffEdge[i]);
353 
354  GetEdgeToElementMap(i,orient[i],mapArray,signArray);
355  for (j=0; j < nmodes[i%2]; j++)
356  {
357  sign = (NekDouble) signArray[j];
358  outarray[ mapArray[j] ] = sign * coeffEdge[i][j];
359  }
360  }
361 
362  int nBoundaryDofs = NumBndryCoeffs();
363  int nInteriorDofs = m_ncoeffs - nBoundaryDofs;
364 
365  if (nInteriorDofs > 0) {
366  Array<OneD, NekDouble> tmp0(m_ncoeffs);
367  Array<OneD, NekDouble> tmp1(m_ncoeffs);
368 
369  StdRegions::StdMatrixKey
370  stdmasskey(StdRegions::eMass,DetShapeType(),*this);
371  MassMatrixOp(outarray,tmp0,stdmasskey);
372  IProductWRTBase(inarray,tmp1);
373 
374  Vmath::Vsub(m_ncoeffs, tmp1, 1, tmp0, 1, tmp1, 1);
375 
376  // get Mass matrix inverse (only of interior DOF)
377  // use block (1,1) of the static condensed system
378  // note: this block alreay contains the inverse matrix
379  MatrixKey
380  masskey(StdRegions::eMass,DetShapeType(),*this);
381  DNekScalMatSharedPtr matsys =
382  (m_staticCondMatrixManager[masskey])->GetBlock(1,1);
383 
384  Array<OneD, NekDouble> rhs(nInteriorDofs);
385  Array<OneD, NekDouble> result(nInteriorDofs);
386 
387  GetInteriorMap(mapArray);
388 
389  for (i = 0; i < nInteriorDofs; i++)
390  {
391  rhs[i] = tmp1[ mapArray[i] ];
392  }
393 
394  Blas::Dgemv('N', nInteriorDofs, nInteriorDofs,
395  matsys->Scale(),
396  &((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 
408 
409  void QuadExp::v_IProductWRTBase(
410  const Array<OneD, const NekDouble>& inarray,
411  Array<OneD, NekDouble> &outarray)
412  {
413  if (m_base[0]->Collocation() && m_base[1]->Collocation())
414  {
415  MultiplyByQuadratureMetric(inarray,outarray);
416  }
417  else
418  {
419  IProductWRTBase_SumFac(inarray,outarray);
420  }
421  }
422 
423 
424  void QuadExp::v_IProductWRTDerivBase(
425  const int dir,
426  const Array<OneD, const NekDouble>& inarray,
427  Array<OneD, NekDouble> & outarray)
428  {
429  IProductWRTDerivBase_SumFac(dir,inarray,outarray);
430  }
431 
432 
433  void QuadExp::v_IProductWRTBase_SumFac(
434  const Array<OneD, const NekDouble>& inarray,
435  Array<OneD, NekDouble> &outarray,
436  bool multiplybyweights)
437  {
438  int nquad0 = m_base[0]->GetNumPoints();
439  int nquad1 = m_base[1]->GetNumPoints();
440  int order0 = m_base[0]->GetNumModes();
441 
442  if(multiplybyweights)
443  {
444  Array<OneD,NekDouble> tmp(nquad0*nquad1+nquad1*order0);
445  Array<OneD,NekDouble> wsp(tmp+nquad0*nquad1);
446 
447  MultiplyByQuadratureMetric(inarray,tmp);
448  StdQuadExp::IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
449  m_base[1]->GetBdata(),
450  tmp,outarray,wsp,true,true);
451  }
452  else
453  {
454  Array<OneD,NekDouble> wsp(nquad1*order0);
455 
456  StdQuadExp::IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
457  m_base[1]->GetBdata(),
458  inarray,outarray,wsp,true,true);
459  }
460  }
461 
462 
463  void QuadExp::v_IProductWRTBase_MatOp(
464  const Array<OneD, const NekDouble>& inarray,
465  Array<OneD, NekDouble> &outarray)
466  {
467  int nq = GetTotPoints();
468  MatrixKey
469  iprodmatkey(StdRegions::eIProductWRTBase,DetShapeType(),*this);
470  DNekScalMatSharedPtr iprodmat = m_matrixManager[iprodmatkey];
471 
472  Blas::Dgemv('N',m_ncoeffs,nq,iprodmat->Scale(),
473  (iprodmat->GetOwnedMatrix())->GetPtr().get(),
474  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
475 
476  }
477 
478 
479  void QuadExp::v_IProductWRTDerivBase_SumFac(
480  const int dir,
481  const Array<OneD, const NekDouble>& inarray,
482  Array<OneD, NekDouble> & outarray)
483  {
484  ASSERTL1((dir==0) || (dir==1) || (dir==2),
485  "Invalid direction.");
486  ASSERTL1((dir==2) ? (m_geom->GetCoordim() ==3):true,
487  "Invalid direction.");
488 
489  int nquad0 = m_base[0]->GetNumPoints();
490  int nquad1 = m_base[1]->GetNumPoints();
491  int nqtot = nquad0*nquad1;
492  int nmodes0 = m_base[0]->GetNumModes();
493 
494  const Array<TwoD, const NekDouble>& df = m_metricinfo->GetDerivFactors(GetPointsKeys());
495 
496  Array<OneD, NekDouble> tmp1(2*nqtot+m_ncoeffs+nmodes0*nquad1);
497  Array<OneD, NekDouble> tmp2(tmp1 + nqtot);
498  Array<OneD, NekDouble> tmp3(tmp1 + 2*nqtot);
499  Array<OneD, NekDouble> tmp4(tmp1 + 2*nqtot+m_ncoeffs);
500 
501  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
502  {
503  Vmath::Vmul(nqtot,
504  &df[2*dir][0], 1,
505  inarray.get(), 1,
506  tmp1.get(), 1);
507  Vmath::Vmul(nqtot,
508  &df[2*dir+1][0], 1,
509  inarray.get(), 1,
510  tmp2.get(),1);
511  }
512  else
513  {
514  Vmath::Smul(nqtot,
515  df[2*dir][0], inarray.get(), 1,
516  tmp1.get(), 1);
517  Vmath::Smul(nqtot,
518  df[2*dir+1][0], inarray.get(), 1,
519  tmp2.get(), 1);
520  }
521 
522  MultiplyByQuadratureMetric(tmp1,tmp1);
523  MultiplyByQuadratureMetric(tmp2,tmp2);
524 
525  IProductWRTBase_SumFacKernel(
526  m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
527  tmp1, tmp3, tmp4, false, true);
528  IProductWRTBase_SumFacKernel(
529  m_base[0]->GetBdata() , m_base[1]->GetDbdata(),
530  tmp2, outarray, tmp4, true, false);
531  Vmath::Vadd(m_ncoeffs, tmp3, 1, outarray, 1, outarray, 1);
532  }
533 
534 
535  void QuadExp::v_IProductWRTDerivBase_MatOp(
536  const int dir,
537  const Array<OneD, const NekDouble>& inarray,
538  Array<OneD, NekDouble> &outarray)
539  {
540  int nq = GetTotPoints();
541  StdRegions::MatrixType mtype = StdRegions::eIProductWRTDerivBase0;
542 
543  switch (dir)
544  {
545  case 0:
546  {
547  mtype = StdRegions::eIProductWRTDerivBase0;
548  }
549  break;
550  case 1:
551  {
552  mtype = StdRegions::eIProductWRTDerivBase1;
553  }
554  break;
555  case 2:
556  {
557  mtype = StdRegions::eIProductWRTDerivBase2;
558  }
559  break;
560  default:
561  {
562  ASSERTL1(false,"input dir is out of range");
563  }
564  break;
565  }
566 
567  MatrixKey iprodmatkey(mtype,DetShapeType(),*this);
568  DNekScalMatSharedPtr iprodmat = m_matrixManager[iprodmatkey];
569 
570  Blas::Dgemv('N', m_ncoeffs, nq, iprodmat->Scale(),
571  (iprodmat->GetOwnedMatrix())->GetPtr().get(),
572  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
573  }
574 
575 
576  void QuadExp::v_NormVectorIProductWRTBase(
577  const Array<OneD, const NekDouble> &Fx,
578  const Array<OneD, const NekDouble> &Fy,
579  const Array<OneD, const NekDouble> &Fz,
580  Array<OneD, NekDouble> &outarray)
581  {
582  int nq = m_base[0]->GetNumPoints()*m_base[1]->GetNumPoints();
583  Array<OneD, NekDouble> Fn(nq);
584 
585  const Array<OneD, const Array<OneD, NekDouble> > &normals =
586  GetLeftAdjacentElementExp()->GetFaceNormal(
587  GetLeftAdjacentElementFace());
588 
589  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
590  {
591  Vmath::Vvtvvtp(nq,&normals[0][0],1,&Fx[0],1,
592  &normals[1][0],1,&Fy[0],1,&Fn[0],1);
593  Vmath::Vvtvp (nq,&normals[2][0],1,&Fz[0],1,&Fn[0],1,&Fn[0],1);
594  }
595  else
596  {
597  Vmath::Svtsvtp(nq,normals[0][0],&Fx[0],1,
598  normals[1][0],&Fy[0],1,&Fn[0],1);
599  Vmath::Svtvp (nq,normals[2][0],&Fz[0],1,&Fn[0],1,&Fn[0],1);
600  }
601 
602  IProductWRTBase(Fn,outarray);
603  }
604 
605  void QuadExp::v_NormVectorIProductWRTBase(
606  const Array<OneD, const Array<OneD, NekDouble> > &Fvec,
607  Array<OneD, NekDouble> &outarray)
608  {
609  NormVectorIProductWRTBase(Fvec[0], Fvec[1], Fvec[2], outarray);
610  }
611 
612  StdRegions::StdExpansionSharedPtr QuadExp::v_GetStdExp(void) const
613  {
614  return MemoryManager<StdRegions::StdQuadExp>
615  ::AllocateSharedPtr(m_base[0]->GetBasisKey(),
616  m_base[1]->GetBasisKey());
617  }
618 
619 
620  StdRegions::StdExpansionSharedPtr QuadExp::v_GetLinStdExp(void) const
621  {
622  LibUtilities::BasisKey bkey0(m_base[0]->GetBasisType(),
623  2, m_base[0]->GetPointsKey());
624  LibUtilities::BasisKey bkey1(m_base[1]->GetBasisType(),
625  2, m_base[1]->GetPointsKey());
626 
627  return MemoryManager<StdRegions::StdQuadExp>
628  ::AllocateSharedPtr( bkey0, bkey1);
629  }
630 
631  void QuadExp::v_GetCoords(
632  Array<OneD, NekDouble> &coords_0,
633  Array<OneD, NekDouble> &coords_1,
634  Array<OneD, NekDouble> &coords_2)
635  {
636  Expansion::v_GetCoords(coords_0, coords_1, coords_2);
637  }
638 
639 
640  void QuadExp::v_GetCoord(const Array<OneD, const NekDouble> &Lcoords,
641  Array<OneD,NekDouble> &coords)
642  {
643  int i;
644 
645  ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[1] <= 1.0 &&
646  Lcoords[1] >= -1.0 && Lcoords[1] <=1.0,
647  "Local coordinates are not in region [-1,1]");
648 
649  m_geom->FillGeom();
650  for (i = 0; i < m_geom->GetCoordim(); ++i)
651  {
652  coords[i] = m_geom->GetCoord(i,Lcoords);
653  }
654  }
655 
656 
657 
658  /**
659  * Given the local cartesian coordinate \a Lcoord evaluate the
660  * value of physvals at this point by calling through to the
661  * StdExpansion method
662  */
663  NekDouble QuadExp::v_StdPhysEvaluate(
664  const Array<OneD, const NekDouble> &Lcoord,
665  const Array<OneD, const NekDouble> &physvals)
666  {
667  // Evaluate point in local (eta) coordinates.
668  return StdQuadExp::v_PhysEvaluate(Lcoord,physvals);
669  }
670 
671  NekDouble QuadExp::v_PhysEvaluate(
672  const Array<OneD, const NekDouble> &coord,
673  const Array<OneD, const NekDouble> &physvals)
674  {
675  Array<OneD,NekDouble> Lcoord = Array<OneD, NekDouble>(2);
676 
677  ASSERTL0(m_geom,"m_geom not defined");
678  m_geom->GetLocCoords(coord,Lcoord);
679 
680  return StdQuadExp::v_PhysEvaluate(Lcoord, physvals);
681  }
682 
683 
684  // Get edge values from the 2D Phys space along an edge
685  // following a counter clockwise edge convention for definition
686  // of edgedir, Note that point distribution is given by QuadExp.
687  void QuadExp::v_GetEdgePhysVals(
688  const int edge,
689  const Array<OneD, const NekDouble> &inarray,
690  Array<OneD,NekDouble> &outarray)
691  {
692  int nquad0 = m_base[0]->GetNumPoints();
693  int nquad1 = m_base[1]->GetNumPoints();
694 
695  StdRegions::Orientation edgedir = GetEorient(edge);
696  switch(edge)
697  {
698  case 0:
699  if (edgedir == StdRegions::eForwards)
700  {
701  Vmath::Vcopy(nquad0,&(inarray[0]),1,&(outarray[0]),1);
702  }
703  else
704  {
705  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0-1),-1,
706  &(outarray[0]),1);
707  }
708  break;
709  case 1:
710  if (edgedir == StdRegions::eForwards)
711  {
712  Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0-1),nquad0,
713  &(outarray[0]),1);
714  }
715  else
716  {
717  Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0*nquad1-1),
718  -nquad0, &(outarray[0]),1);
719  }
720  break;
721  case 2:
722  if (edgedir == StdRegions::eForwards)
723  {
724  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1-1),-1,
725  &(outarray[0]),1);
726  }
727  else
728  {
729  Vmath::Vcopy(nquad0,&(inarray[0])+nquad0*(nquad1-1),1,
730  &(outarray[0]),1);
731  }
732  break;
733  case 3:
734  if (edgedir == StdRegions::eForwards)
735  {
736  Vmath::Vcopy(nquad1,&(inarray[0]) + nquad0*(nquad1-1),
737  -nquad0,&(outarray[0]),1);
738  }
739  else
740  {
741  Vmath::Vcopy(nquad1,&(inarray[0]),nquad0,
742  &(outarray[0]),1);
743  }
744  break;
745  default:
746  ASSERTL0(false,"edge value (< 3) is out of range");
747  break;
748  }
749  }
750 
751 
752  void QuadExp::v_GetTracePhysVals(
753  const int edge,
754  const StdRegions::StdExpansionSharedPtr &EdgeExp,
755  const Array<OneD, const NekDouble> &inarray,
756  Array<OneD,NekDouble> &outarray,
757  StdRegions::Orientation orient)
758  {
759  v_GetEdgePhysVals(edge,EdgeExp,inarray,outarray);
760  }
761 
762 
763  void QuadExp::v_GetEdgePhysVals(
764  const int edge,
765  const StdRegions::StdExpansionSharedPtr &EdgeExp,
766  const Array<OneD, const NekDouble> &inarray,
767  Array<OneD,NekDouble> &outarray)
768  {
769  int nquad0 = m_base[0]->GetNumPoints();
770  int nquad1 = m_base[1]->GetNumPoints();
771 
772  // Implementation for all the basis except Gauss points
773  if (m_base[0]->GetPointsType() !=
774  LibUtilities::eGaussGaussLegendre &&
775  m_base[1]->GetPointsType() !=
776  LibUtilities::eGaussGaussLegendre)
777  {
778  // get points in Cartesian orientation
779  switch (edge)
780  {
781  case 0:
782  Vmath::Vcopy(nquad0,&(inarray[0]),1,&(outarray[0]),1);
783  break;
784  case 1:
785  Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0-1),
786  nquad0,&(outarray[0]),1);
787  break;
788  case 2:
789  Vmath::Vcopy(nquad0,&(inarray[0])+nquad0*(nquad1-1),1,
790  &(outarray[0]),1);
791  break;
792  case 3:
793  Vmath::Vcopy(nquad1,&(inarray[0]),nquad0,
794  &(outarray[0]),1);
795  break;
796  default:
797  ASSERTL0(false,"edge value (< 3) is out of range");
798  break;
799  }
800  }
801  else
802  {
803  QuadExp::v_GetEdgeInterpVals(edge, inarray, outarray);
804  }
805 
806  // Interpolate if required
807  if (m_base[edge%2]->GetPointsKey() !=
808  EdgeExp->GetBasis(0)->GetPointsKey())
809  {
810  Array<OneD,NekDouble> outtmp(max(nquad0,nquad1));
811 
812  outtmp = outarray;
813 
814  LibUtilities::Interp1D(
815  m_base[edge%2]->GetPointsKey(), outtmp,
816  EdgeExp->GetBasis(0)->GetPointsKey(), outarray);
817  }
818 
819  //Reverse data if necessary
820  if(GetCartesianEorient(edge) == StdRegions::eBackwards)
821  {
822  Vmath::Reverse(EdgeExp->GetNumPoints(0),&outarray[0], 1,
823  &outarray[0], 1);
824  }
825  }
826 
827  void QuadExp::v_GetEdgeInterpVals(const int edge,
828  const Array<OneD, const NekDouble> &inarray,
829  Array<OneD,NekDouble> &outarray)
830  {
831  int i;
832  int nq0 = m_base[0]->GetNumPoints();
833  int nq1 = m_base[1]->GetNumPoints();
834 
835  StdRegions::ConstFactorMap factors;
836  factors[StdRegions::eFactorGaussEdge] = edge;
837 
838  StdRegions::StdMatrixKey key(
839  StdRegions::eInterpGauss,
840  DetShapeType(),*this,factors);
841 
842  DNekScalMatSharedPtr mat_gauss = m_matrixManager[key];
843 
844  switch (edge)
845  {
846  case 0:
847  {
848  for (i = 0; i < nq0; i++)
849  {
850  outarray[i] = Blas::Ddot(
851  nq1, mat_gauss->GetOwnedMatrix()->GetPtr().get(),
852  1, &inarray[i], nq0);
853  }
854  break;
855  }
856  case 1:
857  {
858  for (i = 0; i < nq1; i++)
859  {
860  outarray[i] = Blas::Ddot(
861  nq0, mat_gauss->GetOwnedMatrix()->GetPtr().get(),
862  1, &inarray[i * nq0], 1);
863  }
864  break;
865  }
866  case 2:
867  {
868  for (i = 0; i < nq0; i++)
869  {
870  outarray[i] = Blas::Ddot(
871  nq1, mat_gauss->GetOwnedMatrix()->GetPtr().get(),
872  1, &inarray[i], nq0);
873  }
874  break;
875  }
876  case 3:
877  {
878  for (i = 0; i < nq1; i++)
879  {
880  outarray[i] = Blas::Ddot(
881  nq0, mat_gauss->GetOwnedMatrix()->GetPtr().get(),
882  1, &inarray[i * nq0], 1);
883  }
884  break;
885  }
886  default:
887  ASSERTL0(false, "edge value (< 3) is out of range");
888  break;
889  }
890  }
891 
892 
893  void QuadExp::v_GetEdgePhysMap(
894  const int edge,
895  Array<OneD, int> &outarray)
896  {
897  int nquad0 = m_base[0]->GetNumPoints();
898  int nquad1 = m_base[1]->GetNumPoints();
899 
900  // Get points in Cartesian orientation
901  switch (edge)
902  {
903  case 0:
904  outarray = Array<OneD, int>(nquad0);
905  for (int i = 0; i < nquad0; ++i)
906  {
907  outarray[i] = i;
908  }
909  break;
910  case 1:
911  outarray = Array<OneD, int>(nquad1);
912  for (int i = 0; i < nquad1; ++i)
913  {
914  outarray[i] = (nquad0-1) + i*nquad0;
915  }
916  break;
917  case 2:
918  outarray = Array<OneD, int>(nquad0);
919  for (int i = 0; i < nquad0; ++i)
920  {
921  outarray[i] = i + nquad0*(nquad1-1);
922  }
923  break;
924  case 3:
925  outarray = Array<OneD, int>(nquad1);
926  for (int i = 0; i < nquad1; ++i)
927  {
928  outarray[i] = i + i*(nquad0-1);
929  }
930  break;
931  default:
932  ASSERTL0(false, "edge value (< 3) is out of range");
933  break;
934  }
935 
936  }
937 
938 
939 
940 
941  void QuadExp::v_GetEdgeQFactors(
942  const int edge,
943  Array<OneD, NekDouble> &outarray)
944  {
945  int i;
946  int nquad0 = m_base[0]->GetNumPoints();
947  int nquad1 = m_base[1]->GetNumPoints();
948 
949  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
950  const Array<OneD, const NekDouble>& jac = m_metricinfo->GetJac(ptsKeys);
951  const Array<TwoD, const NekDouble>& df = m_metricinfo->GetDerivFactors(ptsKeys);
952 
953  Array<OneD, NekDouble> j (max(nquad0, nquad1), 0.0);
954  Array<OneD, NekDouble> g0(max(nquad0, nquad1), 0.0);
955  Array<OneD, NekDouble> g1(max(nquad0, nquad1), 0.0);
956  Array<OneD, NekDouble> g2(max(nquad0, nquad1), 0.0);
957  Array<OneD, NekDouble> g3(max(nquad0, nquad1), 0.0);
958 
959  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
960  {
961  // Implementation for all the basis except Gauss points
962  if (m_base[0]->GetPointsType()
963  != LibUtilities::eGaussGaussLegendre
964  && m_base[1]->GetPointsType() !=
965  LibUtilities::eGaussGaussLegendre)
966  {
967  switch (edge)
968  {
969  case 0:
970  Vmath::Vcopy(nquad0, &(df[1][0]),
971  1, &(g1[0]), 1);
972  Vmath::Vcopy(nquad0, &(df[3][0]),
973  1, &(g3[0]), 1);
974  Vmath::Vcopy(nquad0, &(jac[0]),1, &(j[0]), 1);
975 
976  for (i = 0; i < nquad0; ++i)
977  {
978  outarray[i] = j[i]*sqrt(g1[i]*g1[i]
979  + g3[i]*g3[i]);
980  }
981  break;
982  case 1:
983  Vmath::Vcopy(nquad1,
984  &(df[0][0])+(nquad0-1), nquad0,
985  &(g0[0]), 1);
986 
987  Vmath::Vcopy(nquad1,
988  &(df[2][0])+(nquad0-1), nquad0,
989  &(g2[0]), 1);
990 
991  Vmath::Vcopy(nquad1,
992  &(jac[0])+(nquad0-1), nquad0,
993  &(j[0]), 1);
994 
995  for (i = 0; i < nquad1; ++i)
996  {
997  outarray[i] = j[i]*sqrt(g0[i]*g0[i] +
998  g2[i]*g2[i]);
999  }
1000  break;
1001  case 2:
1002 
1003  Vmath::Vcopy(nquad0,
1004  &(df[1][0])+(nquad0*nquad1-1), -1,
1005  &(g1[0]), 1);
1006 
1007  Vmath::Vcopy(nquad0,
1008  &(df[3][0])+(nquad0*nquad1-1), -1,
1009  &(g3[0]), 1);
1010 
1011  Vmath::Vcopy(nquad0,
1012  &(jac[0])+(nquad0*nquad1-1), -1,
1013  &(j[0]), 1);
1014 
1015  for (i = 0; i < nquad0; ++i)
1016  {
1017  outarray[i] = j[i]*sqrt(g1[i]*g1[i]
1018  + g3[i]*g3[i]);
1019  }
1020  break;
1021  case 3:
1022 
1023  Vmath::Vcopy(nquad1,
1024  &(df[0][0])+nquad0*(nquad1-1),
1025  -nquad0,&(g0[0]), 1);
1026 
1027  Vmath::Vcopy(nquad1,
1028  &(df[2][0])+nquad0*(nquad1-1),
1029  -nquad0,&(g2[0]), 1);
1030 
1031  Vmath::Vcopy(nquad1,
1032  &(jac[0])+nquad0*(nquad1-1), -nquad0,
1033  &(j[0]), 1);
1034 
1035  for (i = 0; i < nquad1; ++i)
1036  {
1037  outarray[i] = j[i]*sqrt(g0[i]*g0[i] +
1038  g2[i]*g2[i]);
1039  }
1040  break;
1041  default:
1042  ASSERTL0(false,"edge value (< 3) is out of range");
1043  break;
1044  }
1045  }
1046  else
1047  {
1048  int nqtot = nquad0 * nquad1;
1049  Array<OneD, NekDouble> tmp_gmat0(nqtot, 0.0);
1050  Array<OneD, NekDouble> tmp_gmat1(nqtot, 0.0);
1051  Array<OneD, NekDouble> tmp_gmat2(nqtot, 0.0);
1052  Array<OneD, NekDouble> tmp_gmat3(nqtot, 0.0);
1053  Array<OneD, NekDouble> g0_edge(max(nquad0, nquad1), 0.0);
1054  Array<OneD, NekDouble> g1_edge(max(nquad0, nquad1), 0.0);
1055  Array<OneD, NekDouble> g2_edge(max(nquad0, nquad1), 0.0);
1056  Array<OneD, NekDouble> g3_edge(max(nquad0, nquad1), 0.0);
1057  Array<OneD, NekDouble> jac_edge(max(nquad0, nquad1), 0.0);
1058 
1059  switch (edge)
1060  {
1061  case 0:
1062  Vmath::Vmul(nqtot,&(df[1][0]),1,&jac[0],1,
1063  &(tmp_gmat1[0]),1);
1064  Vmath::Vmul(nqtot,&(df[3][0]),1,&jac[0],1,
1065  &(tmp_gmat3[0]),1);
1066  QuadExp::v_GetEdgeInterpVals(
1067  edge, tmp_gmat1, g1_edge);
1068  QuadExp::v_GetEdgeInterpVals(
1069  edge, tmp_gmat3, g3_edge);
1070 
1071  for (i = 0; i < nquad0; ++i)
1072  {
1073  outarray[i] = sqrt(g1_edge[i]*g1_edge[i] +
1074  g3_edge[i]*g3_edge[i]);
1075  }
1076  break;
1077 
1078  case 1:
1079  Vmath::Vmul(nqtot,
1080  &(df[0][0]), 1,
1081  &jac[0], 1,
1082  &(tmp_gmat0[0]), 1);
1083  Vmath::Vmul(nqtot,
1084  &(df[2][0]), 1,
1085  &jac[0], 1,
1086  &(tmp_gmat2[0]),
1087  1);
1088  QuadExp::v_GetEdgeInterpVals(
1089  edge, tmp_gmat0, g0_edge);
1090  QuadExp::v_GetEdgeInterpVals(
1091  edge, tmp_gmat2, g2_edge);
1092 
1093  for (i = 0; i < nquad1; ++i)
1094  {
1095  outarray[i] = sqrt(g0_edge[i]*g0_edge[i]
1096  + g2_edge[i]*g2_edge[i]);
1097  }
1098 
1099  break;
1100  case 2:
1101 
1102  Vmath::Vmul(nqtot,
1103  &(df[1][0]), 1,
1104  &jac[0], 1,
1105  &(tmp_gmat1[0]), 1);
1106  Vmath::Vmul(nqtot,
1107  &(df[3][0]), 1,
1108  &jac[0], 1,
1109  &(tmp_gmat3[0]),1);
1110  QuadExp::v_GetEdgeInterpVals(
1111  edge, tmp_gmat1, g1_edge);
1112  QuadExp::v_GetEdgeInterpVals(
1113  edge, tmp_gmat3, g3_edge);
1114 
1115 
1116  for (i = 0; i < nquad0; ++i)
1117  {
1118  outarray[i] = sqrt(g1_edge[i]*g1_edge[i]
1119  + g3_edge[i]*g3_edge[i]);
1120  }
1121 
1122  Vmath::Reverse(nquad0,&outarray[0],1,&outarray[0],1);
1123 
1124  break;
1125  case 3:
1126  Vmath::Vmul(nqtot,
1127  &(df[0][0]), 1,
1128  &jac[0], 1,
1129  &(tmp_gmat0[0]), 1);
1130  Vmath::Vmul(nqtot,
1131  &(df[2][0]),1,
1132  &jac[0], 1,
1133  &(tmp_gmat2[0]),1);
1134  QuadExp::v_GetEdgeInterpVals(
1135  edge, tmp_gmat0, g0_edge);
1136  QuadExp::v_GetEdgeInterpVals(
1137  edge, tmp_gmat2, g2_edge);
1138 
1139 
1140  for (i = 0; i < nquad1; ++i)
1141  {
1142  outarray[i] = sqrt(g0_edge[i]*g0_edge[i] +
1143  g2_edge[i]*g2_edge[i]);
1144  }
1145 
1146  Vmath::Reverse(nquad1,
1147  &outarray[0], 1,
1148  &outarray[0], 1);
1149 
1150  break;
1151  default:
1152  ASSERTL0(false,"edge value (< 3) is out of range");
1153  break;
1154  }
1155  }
1156  }
1157  else
1158  {
1159 
1160  switch (edge)
1161  {
1162  case 0:
1163 
1164 
1165 
1166  for (i = 0; i < nquad0; ++i)
1167  {
1168  outarray[i] = jac[0]*sqrt(df[1][0]*df[1][0] +
1169  df[3][0]*df[3][0]);
1170  }
1171  break;
1172  case 1:
1173  for (i = 0; i < nquad1; ++i)
1174  {
1175  outarray[i] = jac[0]*sqrt(df[0][0]*df[0][0] +
1176  df[2][0]*df[2][0]);
1177  }
1178  break;
1179  case 2:
1180  for (i = 0; i < nquad0; ++i)
1181  {
1182  outarray[i] = jac[0]*sqrt(df[1][0]*df[1][0] +
1183  df[3][0]*df[3][0]);
1184  }
1185  break;
1186  case 3:
1187  for (i = 0; i < nquad1; ++i)
1188  {
1189  outarray[i] = jac[0]*sqrt(df[0][0]*df[0][0] +
1190  df[2][0]*df[2][0]);
1191  }
1192  break;
1193  default:
1194  ASSERTL0(false,"edge value (< 3) is out of range");
1195  break;
1196  }
1197  }
1198  }
1199 
1200 
1201  void QuadExp::v_ComputeEdgeNormal(const int edge)
1202  {
1203  int i;
1204  const SpatialDomains::GeomFactorsSharedPtr & geomFactors =
1205  GetGeom()->GetMetricInfo();
1206  SpatialDomains::GeomType type = geomFactors->GetGtype();
1207  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
1208  const Array<TwoD, const NekDouble> & df = geomFactors->GetDerivFactors(ptsKeys);
1209  const Array<OneD, const NekDouble> & jac = geomFactors->GetJac(ptsKeys);
1210  int nqe;
1211  if (edge == 0 || edge == 2)
1212  {
1213  nqe = m_base[0]->GetNumPoints();
1214  }
1215  else
1216  {
1217  nqe = m_base[1]->GetNumPoints();
1218  }
1219  int vCoordDim = GetCoordim();
1220 
1221  m_edgeNormals[edge] = Array<OneD, Array<OneD, NekDouble> >
1222  (vCoordDim);
1223  Array<OneD, Array<OneD, NekDouble> > &normal = m_edgeNormals[edge];
1224  for (i = 0; i < vCoordDim; ++i)
1225  {
1226  normal[i] = Array<OneD, NekDouble>(nqe);
1227  }
1228 
1229  // Regular geometry case
1230  if ((type == SpatialDomains::eRegular)||
1231  (type == SpatialDomains::eMovingRegular))
1232  {
1233  NekDouble fac;
1234  // Set up normals
1235  switch (edge)
1236  {
1237  case 0:
1238  for (i = 0; i < vCoordDim; ++i)
1239  {
1240  Vmath::Fill(nqe, -df[2*i+1][0], normal[i], 1);
1241  }
1242  break;
1243  case 1:
1244  for (i = 0; i < vCoordDim; ++i)
1245  {
1246  Vmath::Fill(nqe, df[2*i][0], normal[i], 1);
1247  }
1248  break;
1249  case 2:
1250  for (i = 0; i < vCoordDim; ++i)
1251  {
1252  Vmath::Fill(nqe, df[2*i+1][0], normal[i], 1);
1253  }
1254  break;
1255  case 3:
1256  for (i = 0; i < vCoordDim; ++i)
1257  {
1258  Vmath::Fill(nqe, -df[2*i][0], normal[i], 1);
1259  }
1260  break;
1261  default:
1262  ASSERTL0(false, "edge is out of range (edge < 4)");
1263  }
1264 
1265  // normalise
1266  fac = 0.0;
1267  for (i =0 ; i < vCoordDim; ++i)
1268  {
1269  fac += normal[i][0]*normal[i][0];
1270  }
1271  fac = 1.0/sqrt(fac);
1272  for (i = 0; i < vCoordDim; ++i)
1273  {
1274  Vmath::Smul(nqe, fac, normal[i], 1,normal[i], 1);
1275  }
1276  }
1277  else // Set up deformed normals
1278  {
1279  int j;
1280 
1281  int nquad0 = ptsKeys[0].GetNumPoints();
1282  int nquad1 = ptsKeys[1].GetNumPoints();
1283 
1284  LibUtilities::PointsKey from_key;
1285 
1286  Array<OneD,NekDouble> normals(vCoordDim*max(nquad0,nquad1),0.0);
1287  Array<OneD,NekDouble> edgejac(vCoordDim*max(nquad0,nquad1),0.0);
1288 
1289  // Extract Jacobian along edges and recover local
1290  // derivates (dx/dr) for polynomial interpolation by
1291  // multiplying m_gmat by jacobian
1292 
1293  // Implementation for all the basis except Gauss points
1294  if (m_base[0]->GetPointsType() !=
1295  LibUtilities::eGaussGaussLegendre
1296  && m_base[1]->GetPointsType() !=
1297  LibUtilities::eGaussGaussLegendre)
1298  {
1299  switch (edge)
1300  {
1301  case 0:
1302  for (j = 0; j < nquad0; ++j)
1303  {
1304  edgejac[j] = jac[j];
1305  for (i = 0; i < vCoordDim; ++i)
1306  {
1307  normals[i*nquad0+j] =
1308  -df[2*i+1][j]*edgejac[j];
1309  }
1310  }
1311  from_key = ptsKeys[0];
1312  break;
1313  case 1:
1314  for (j = 0; j < nquad1; ++j)
1315  {
1316  edgejac[j] = jac[nquad0*j+nquad0-1];
1317  for (i = 0; i < vCoordDim; ++i)
1318  {
1319  normals[i*nquad1+j] =
1320  df[2*i][nquad0*j + nquad0-1]
1321  *edgejac[j];
1322  }
1323  }
1324  from_key = ptsKeys[1];
1325  break;
1326  case 2:
1327  for (j = 0; j < nquad0; ++j)
1328  {
1329  edgejac[j] = jac[nquad0*(nquad1-1)+j];
1330  for (i = 0; i < vCoordDim; ++i)
1331  {
1332  normals[i*nquad0+j] =
1333  (df[2*i+1][nquad0*(nquad1-1)+j])
1334  *edgejac[j];
1335  }
1336  }
1337  from_key = ptsKeys[0];
1338  break;
1339  case 3:
1340  for (j = 0; j < nquad1; ++j)
1341  {
1342  edgejac[j] = jac[nquad0*j];
1343  for (i = 0; i < vCoordDim; ++i)
1344  {
1345  normals[i*nquad1+j] =
1346  -df[2*i][nquad0*j]*edgejac[j];
1347  }
1348  }
1349  from_key = ptsKeys[1];
1350  break;
1351  default:
1352  ASSERTL0(false,"edge is out of range (edge < 3)");
1353  }
1354  }
1355  else
1356  {
1357  int nqtot = nquad0 * nquad1;
1358  Array<OneD, NekDouble> tmp_gmat(nqtot, 0.0);
1359  Array<OneD, NekDouble> tmp_gmat_edge(nqe, 0.0);
1360 
1361  switch (edge)
1362  {
1363  case 0:
1364  for (j = 0; j < nquad0; ++j)
1365  {
1366  for (i = 0; i < vCoordDim; ++i)
1367  {
1368  Vmath::Vmul(nqtot,
1369  &(df[2*i+1][0]), 1,
1370  &jac[0], 1,
1371  &(tmp_gmat[0]), 1);
1372  QuadExp::v_GetEdgeInterpVals(
1373  edge, tmp_gmat, tmp_gmat_edge);
1374  normals[i*nquad0+j] = -tmp_gmat_edge[j];
1375  }
1376  }
1377  from_key = ptsKeys[0];
1378  break;
1379  case 1:
1380  for (j = 0; j < nquad1; ++j)
1381  {
1382  for (i = 0; i < vCoordDim; ++i)
1383  {
1384  Vmath::Vmul(nqtot,
1385  &(df[2*i][0]), 1,
1386  &jac[0], 1,
1387  &(tmp_gmat[0]), 1);
1388  QuadExp::v_GetEdgeInterpVals(
1389  edge, tmp_gmat, tmp_gmat_edge);
1390  normals[i*nquad1+j] = tmp_gmat_edge[j];
1391  }
1392  }
1393  from_key = ptsKeys[1];
1394  break;
1395  case 2:
1396  for (j = 0; j < nquad0; ++j)
1397  {
1398  for (i = 0; i < vCoordDim; ++i)
1399  {
1400  Vmath::Vmul(nqtot,
1401  &(df[2*i+1][0]), 1,
1402  &jac[0], 1,
1403  &(tmp_gmat[0]), 1);
1404  QuadExp::v_GetEdgeInterpVals(
1405  edge, tmp_gmat, tmp_gmat_edge);
1406  normals[i*nquad0+j] = tmp_gmat_edge[j];
1407  }
1408  }
1409  from_key = ptsKeys[0];
1410  break;
1411  case 3:
1412  for (j = 0; j < nquad1; ++j)
1413  {
1414  for (i = 0; i < vCoordDim; ++i)
1415  {
1416  Vmath::Vmul(nqtot,
1417  &(df[2*i][0]), 1,
1418  &jac[0], 1,
1419  &(tmp_gmat[0]) ,1);
1420  QuadExp::v_GetEdgeInterpVals(
1421  edge, tmp_gmat, tmp_gmat_edge);
1422  normals[i*nquad1+j] = -tmp_gmat_edge[j];
1423  }
1424  }
1425  from_key = ptsKeys[1];
1426  break;
1427  default:
1428  ASSERTL0(false,"edge is out of range (edge < 3)");
1429  }
1430  }
1431 
1432  int nq = from_key.GetNumPoints();
1433  Array<OneD,NekDouble> work(nqe,0.0);
1434 
1435  // interpolate Jacobian and invert
1436  LibUtilities::Interp1D(
1437  from_key,jac, m_base[0]->GetPointsKey(), work);
1438  Vmath::Sdiv(nq,1.0,&work[0],1,&work[0],1);
1439 
1440  // interpolate
1441  for (i = 0; i < GetCoordim(); ++i)
1442  {
1443  LibUtilities::Interp1D(
1444  from_key,&normals[i*nq],
1445  m_base[0]->GetPointsKey(),
1446  &normal[i][0]);
1447  Vmath::Vmul(nqe, work, 1, normal[i], 1, normal[i], 1);
1448  }
1449 
1450  //normalise normal vectors
1451  Vmath::Zero(nqe,work,1);
1452  for (i = 0; i < GetCoordim(); ++i)
1453  {
1454  Vmath::Vvtvp(nqe,
1455  normal[i], 1,
1456  normal[i],1 ,
1457  work, 1,
1458  work, 1);
1459  }
1460 
1461  Vmath::Vsqrt(nqe,work,1,work,1);
1462  Vmath::Sdiv(nqe,1.0,work,1,work,1);
1463 
1464  for (i = 0; i < GetCoordim(); ++i)
1465  {
1466  Vmath::Vmul(nqe, normal[i], 1, work, 1, normal[i], 1);
1467  }
1468 
1469  // Reverse direction so that points are in
1470  // anticlockwise direction if edge >=2
1471  if (edge >= 2)
1472  {
1473  for (i = 0; i < GetCoordim(); ++i)
1474  {
1475  Vmath::Reverse(nqe, normal[i], 1, normal[i], 1);
1476  }
1477  }
1478  }
1479  if (GetGeom()->GetEorient(edge) == StdRegions::eBackwards)
1480  {
1481  for (i = 0; i < vCoordDim; ++i)
1482  {
1483  if (geomFactors->GetGtype() == SpatialDomains::eDeformed)
1484  {
1485  Vmath::Reverse(nqe, normal[i], 1, normal[i],1);
1486  }
1487  }
1488  }
1489  }
1490 
1491  const SpatialDomains::GeomFactorsSharedPtr& QuadExp::v_GetMetricInfo() const
1492  {
1493  return m_metricinfo;
1494  }
1495 
1496 
1497  int QuadExp::v_GetCoordim()
1498  {
1499  return m_geom->GetCoordim();
1500  }
1501 
1502 
1503  void QuadExp::v_ExtractDataToCoeffs(
1504  const NekDouble *data,
1505  const std::vector<unsigned int > &nummodes,
1506  int mode_offset,
1507  NekDouble *coeffs,
1508  std::vector<LibUtilities::BasisType> &fromType)
1509  {
1510  int data_order0 = nummodes[mode_offset];
1511  int fillorder0 = std::min(m_base[0]->GetNumModes(),data_order0);
1512 
1513  int data_order1 = nummodes[mode_offset + 1];
1514  int order1 = m_base[1]->GetNumModes();
1515  int fillorder1 = min(order1,data_order1);
1516 
1517  // Check if same basis
1518  if (fromType[0] != m_base[0]->GetBasisType() ||
1519  fromType[1] != m_base[1]->GetBasisType())
1520  {
1521  // Construct a quad with the appropriate basis type at our
1522  // quadrature points, and one more to do a forwards
1523  // transform. We can then copy the output to coeffs.
1524  StdRegions::StdQuadExp tmpQuad(
1525  LibUtilities::BasisKey(
1526  fromType[0], data_order0, m_base[0]->GetPointsKey()),
1527  LibUtilities::BasisKey(
1528  fromType[1], data_order1, m_base[1]->GetPointsKey()));
1529  StdRegions::StdQuadExp tmpQuad2(m_base[0]->GetBasisKey(),
1530  m_base[1]->GetBasisKey());
1531 
1532  Array<OneD, const NekDouble> tmpData(tmpQuad.GetNcoeffs(), data);
1533  Array<OneD, NekDouble> tmpBwd(tmpQuad2.GetTotPoints());
1534  Array<OneD, NekDouble> tmpOut(tmpQuad2.GetNcoeffs());
1535 
1536  tmpQuad.BwdTrans(tmpData, tmpBwd);
1537  tmpQuad2.FwdTrans(tmpBwd, tmpOut);
1538  Vmath::Vcopy(tmpOut.num_elements(), &tmpOut[0], 1, coeffs, 1);
1539 
1540  return;
1541  }
1542 
1543  switch (m_base[0]->GetBasisType())
1544  {
1545  case LibUtilities::eModified_A:
1546  {
1547  int i;
1548  int cnt = 0;
1549  int cnt1 = 0;
1550 
1551  ASSERTL1(m_base[1]->GetBasisType() ==
1552  LibUtilities::eModified_A,
1553  "Extraction routine not set up for this basis");
1554 
1555  Vmath::Zero(m_ncoeffs,coeffs,1);
1556  for (i = 0; i < fillorder0; ++i)
1557  {
1558  Vmath::Vcopy(fillorder1, data + cnt, 1, coeffs +cnt1, 1);
1559  cnt += data_order1;
1560  cnt1 += order1;
1561  }
1562  }
1563  break;
1564  case LibUtilities::eGLL_Lagrange:
1565  {
1566  LibUtilities::PointsKey
1567  p0(nummodes[0], LibUtilities::eGaussLobattoLegendre);
1568  LibUtilities::PointsKey
1569  p1(nummodes[1], LibUtilities::eGaussLobattoLegendre);
1570  LibUtilities::PointsKey t0(
1571  m_base[0]->GetNumModes(),
1572  LibUtilities::eGaussLobattoLegendre);
1573  LibUtilities::PointsKey t1(
1574  m_base[1]->GetNumModes(),
1575  LibUtilities::eGaussLobattoLegendre);
1576  LibUtilities::Interp2D(p0, p1, data, t0, t1, coeffs);
1577  }
1578  break;
1579  case LibUtilities::eGauss_Lagrange:
1580  {
1581  // Assume that input is also Gll_Lagrange but no way to check;
1582  LibUtilities::PointsKey
1583  p0(nummodes[0],LibUtilities::eGaussGaussLegendre);
1584  LibUtilities::PointsKey
1585  p1(nummodes[1],LibUtilities::eGaussGaussLegendre);
1586  LibUtilities::PointsKey t0(
1587  m_base[0]->GetNumModes(),
1588  LibUtilities::eGaussGaussLegendre);
1589  LibUtilities::PointsKey t1(
1590  m_base[1]->GetNumModes(),
1591  LibUtilities::eGaussGaussLegendre);
1592  LibUtilities::Interp2D(p0, p1, data, t0, t1, coeffs);
1593  }
1594  break;
1595  default:
1596  ASSERTL0(false,
1597  "basis is either not set up or not hierarchicial");
1598  }
1599  }
1600 
1601 
1602  StdRegions::Orientation QuadExp::v_GetEorient(int edge)
1603  {
1604  return m_geom->GetEorient(edge);
1605  }
1606 
1607 
1608  StdRegions::Orientation QuadExp::v_GetCartesianEorient(int edge)
1609  {
1610  return GetGeom2D()->GetCartesianEorient(edge);
1611  }
1612 
1613 
1614  const LibUtilities::BasisSharedPtr& QuadExp::v_GetBasis(int dir) const
1615  {
1616  ASSERTL1(dir >= 0 &&dir <= 1, "input dir is out of range");
1617  return m_base[dir];
1618  }
1619 
1620 
1621  int QuadExp::v_GetNumPoints(const int dir) const
1622  {
1623  return GetNumPoints(dir);
1624  }
1625 
1626  DNekMatSharedPtr QuadExp::v_GenMatrix(
1627  const StdRegions::StdMatrixKey &mkey)
1628  {
1629  DNekMatSharedPtr returnval;
1630  switch (mkey.GetMatrixType())
1631  {
1632  case StdRegions::eHybridDGHelmholtz:
1633  case StdRegions::eHybridDGLamToU:
1634  case StdRegions::eHybridDGLamToQ0:
1635  case StdRegions::eHybridDGLamToQ1:
1636  case StdRegions::eHybridDGLamToQ2:
1637  case StdRegions::eHybridDGHelmBndLam:
1638  case StdRegions::eInvLaplacianWithUnityMean:
1639  returnval = Expansion2D::v_GenMatrix(mkey);
1640  break;
1641  default:
1642  returnval = StdQuadExp::v_GenMatrix(mkey);
1643  }
1644  return returnval;
1645  }
1646 
1647  DNekMatSharedPtr QuadExp::v_CreateStdMatrix(
1648  const StdRegions::StdMatrixKey &mkey)
1649  {
1650  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
1651  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
1652  StdRegions::StdQuadExpSharedPtr tmp =
1653  MemoryManager<StdQuadExp>::AllocateSharedPtr(bkey0,bkey1);
1654  return tmp->GetStdMatrix(mkey);
1655  }
1656 
1657 
1658  DNekScalMatSharedPtr QuadExp::CreateMatrix(const MatrixKey &mkey)
1659  {
1660  DNekScalMatSharedPtr returnval;
1661  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
1662 
1663  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,
1664  "Geometric information is not set up");
1665 
1666  switch (mkey.GetMatrixType())
1667  {
1668  case StdRegions::eMass:
1669  {
1670  if ((m_metricinfo->GetGtype() ==
1671  SpatialDomains::eDeformed) || (mkey.GetNVarCoeff()))
1672  {
1673  NekDouble one = 1.0;
1674  DNekMatSharedPtr mat = GenMatrix(mkey);
1675  returnval = MemoryManager<DNekScalMat>::
1676  AllocateSharedPtr(one,mat);
1677  }
1678  else
1679  {
1680  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1681  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1682  returnval = MemoryManager<DNekScalMat>::
1683  AllocateSharedPtr(jac,mat);
1684  }
1685  }
1686  break;
1687  case StdRegions::eInvMass:
1688  {
1689  if ((m_metricinfo->GetGtype() ==
1690  SpatialDomains::eDeformed) || (mkey.GetNVarCoeff()))
1691  {
1692  NekDouble one = 1.0;
1693  StdRegions::StdMatrixKey masskey(
1694  StdRegions::eMass, DetShapeType(), *this);
1695  DNekMatSharedPtr mat = GenMatrix(masskey);
1696  mat->Invert();
1697 
1698  returnval = MemoryManager<DNekScalMat>::
1699  AllocateSharedPtr(one,mat);
1700  }
1701  else
1702  {
1703  NekDouble fac = 1.0/(m_metricinfo->GetJac(ptsKeys))[0];
1704  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1705  returnval = MemoryManager<DNekScalMat>::
1706  AllocateSharedPtr(fac,mat);
1707  }
1708  }
1709  break;
1710  case StdRegions::eWeakDeriv0:
1711  case StdRegions::eWeakDeriv1:
1712  case StdRegions::eWeakDeriv2:
1713  {
1714  if((m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1715  || (mkey.GetNVarCoeff()))
1716  {
1717  NekDouble one = 1.0;
1718  DNekMatSharedPtr mat = GenMatrix(mkey);
1719 
1720  returnval = MemoryManager<DNekScalMat>::
1721  AllocateSharedPtr(one,mat);
1722  }
1723  else
1724  {
1725  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1726  Array<TwoD, const NekDouble> df =
1727  m_metricinfo->GetDerivFactors(ptsKeys);
1728  int dir = 0;
1729 
1730  switch(mkey.GetMatrixType())
1731  {
1732  case StdRegions::eWeakDeriv0:
1733  dir = 0;
1734  break;
1735  case StdRegions::eWeakDeriv1:
1736  dir = 1;
1737  break;
1738  case StdRegions::eWeakDeriv2:
1739  dir = 2;
1740  break;
1741  default:
1742  break;
1743  }
1744 
1745  MatrixKey deriv0key(StdRegions::eWeakDeriv0,
1746  mkey.GetShapeType(), *this);
1747  MatrixKey deriv1key(StdRegions::eWeakDeriv1,
1748  mkey.GetShapeType(), *this);
1749 
1750  DNekMat &deriv0 = *GetStdMatrix(deriv0key);
1751  DNekMat &deriv1 = *GetStdMatrix(deriv1key);
1752 
1753  int rows = deriv0.GetRows();
1754  int cols = deriv1.GetColumns();
1755 
1756  DNekMatSharedPtr WeakDeriv = MemoryManager<DNekMat>::
1757  AllocateSharedPtr(rows,cols);
1758  (*WeakDeriv) = df[2*dir][0]*deriv0 +
1759  df[2*dir+1][0]*deriv1;
1760  returnval = MemoryManager<DNekScalMat>::
1761  AllocateSharedPtr(jac,WeakDeriv);
1762  }
1763  }
1764  break;
1765  case StdRegions::eLaplacian:
1766  {
1767  if( (m_metricinfo->GetGtype() ==
1768  SpatialDomains::eDeformed) || (mkey.GetNVarCoeff() > 0)
1769  || (mkey.ConstFactorExists
1770  (StdRegions::eFactorSVVCutoffRatio)))
1771  {
1772  NekDouble one = 1.0;
1773  DNekMatSharedPtr mat = GenMatrix(mkey);
1774 
1775  returnval = MemoryManager<DNekScalMat>::
1776  AllocateSharedPtr(one,mat);
1777  }
1778  else
1779  {
1780  MatrixKey lap00key(StdRegions::eLaplacian00,
1781  mkey.GetShapeType(), *this);
1782  MatrixKey lap01key(StdRegions::eLaplacian01,
1783  mkey.GetShapeType(), *this);
1784  MatrixKey lap11key(StdRegions::eLaplacian11,
1785  mkey.GetShapeType(), *this);
1786 
1787  DNekMat &lap00 = *GetStdMatrix(lap00key);
1788  DNekMat &lap01 = *GetStdMatrix(lap01key);
1789  DNekMat &lap11 = *GetStdMatrix(lap11key);
1790 
1791  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1792  Array<TwoD, const NekDouble>
1793  gmat = m_metricinfo->GetGmat(ptsKeys);
1794 
1795  int rows = lap00.GetRows();
1796  int cols = lap00.GetColumns();
1797 
1798  DNekMatSharedPtr lap =
1799  MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols);
1800 
1801  (*lap) = gmat[0][0] * lap00 +
1802  gmat[1][0] * (lap01 + Transpose(lap01)) +
1803  gmat[3][0] * lap11;
1804 
1805  returnval = MemoryManager<DNekScalMat>::
1806  AllocateSharedPtr(jac,lap);
1807  }
1808  }
1809  break;
1810  case StdRegions::eInvLaplacianWithUnityMean:
1811  {
1812  DNekMatSharedPtr mat = GenMatrix(mkey);
1813  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0,mat);
1814  }
1815  break;
1816  case StdRegions::eHelmholtz:
1817  {
1818  NekDouble lambda =
1819  mkey.GetConstFactor(StdRegions::eFactorLambda);
1820 
1821  MatrixKey masskey(mkey, StdRegions::eMass);
1822  DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
1823 
1824  MatrixKey lapkey(mkey, StdRegions::eLaplacian);
1825  DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
1826 
1827  int rows = LapMat.GetRows();
1828  int cols = LapMat.GetColumns();
1829 
1830  DNekMatSharedPtr helm = MemoryManager<DNekMat>::
1831  AllocateSharedPtr(rows,cols);
1832 
1833  NekDouble one = 1.0;
1834  (*helm) = LapMat + lambda*MassMat;
1835 
1836  returnval =
1837  MemoryManager<DNekScalMat>::AllocateSharedPtr(one,helm);
1838  }
1839  break;
1840  case StdRegions::eIProductWRTBase:
1841  {
1842  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1843  {
1844  NekDouble one = 1.0;
1845  DNekMatSharedPtr mat = GenMatrix(mkey);
1846  returnval = MemoryManager<DNekScalMat>::
1847  AllocateSharedPtr(one,mat);
1848  }
1849  else
1850  {
1851  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1852  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1853  returnval = MemoryManager<DNekScalMat>::
1854  AllocateSharedPtr(jac,mat);
1855  }
1856  }
1857  break;
1858  case StdRegions::eIProductWRTDerivBase0:
1859  case StdRegions::eIProductWRTDerivBase1:
1860  case StdRegions::eIProductWRTDerivBase2:
1861  {
1862  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1863  {
1864  NekDouble one = 1.0;
1865  DNekMatSharedPtr mat = GenMatrix(mkey);
1866  returnval = MemoryManager<DNekScalMat>::
1867  AllocateSharedPtr(one,mat);
1868  }
1869  else
1870  {
1871  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1872  const Array<TwoD, const NekDouble>& df =
1873  m_metricinfo->GetDerivFactors(ptsKeys);
1874  int dir = 0;
1875 
1876  switch(mkey.GetMatrixType())
1877  {
1878  case StdRegions::eIProductWRTDerivBase0:
1879  dir = 0;
1880  break;
1881  case StdRegions::eIProductWRTDerivBase1:
1882  dir = 1;
1883  break;
1884  case StdRegions::eIProductWRTDerivBase2:
1885  dir = 2;
1886  break;
1887  default:
1888  break;
1889  }
1890 
1891  MatrixKey iProdDeriv0Key(
1892  StdRegions::eIProductWRTDerivBase0,
1893  mkey.GetShapeType(), *this);
1894  MatrixKey iProdDeriv1Key(
1895  StdRegions::eIProductWRTDerivBase1,
1896  mkey.GetShapeType(), *this);
1897 
1898  DNekMat &stdiprod0 = *GetStdMatrix(iProdDeriv0Key);
1899  DNekMat &stdiprod1 = *GetStdMatrix(iProdDeriv0Key);
1900 
1901  int rows = stdiprod0.GetRows();
1902  int cols = stdiprod1.GetColumns();
1903 
1904  DNekMatSharedPtr mat = MemoryManager<DNekMat>::
1905  AllocateSharedPtr(rows,cols);
1906  (*mat) = df[2*dir][0]*stdiprod0 +
1907  df[2*dir+1][0]*stdiprod1;
1908 
1909  returnval = MemoryManager<DNekScalMat>::
1910  AllocateSharedPtr(jac,mat);
1911  }
1912  }
1913  break;
1914  case StdRegions::eInvHybridDGHelmholtz:
1915  {
1916  NekDouble one = 1.0;
1917 
1918  MatrixKey hkey(StdRegions::eHybridDGHelmholtz,
1919  DetShapeType(), *this,
1920  mkey.GetConstFactors(), mkey.GetVarCoeffs());
1921  DNekMatSharedPtr mat = GenMatrix(hkey);
1922 
1923  mat->Invert();
1924  returnval =
1925  MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1926  }
1927  break;
1928  case StdRegions::eInterpGauss:
1929  {
1930  DNekMatSharedPtr m_Ix;
1931  Array<OneD, NekDouble> coords(1, 0.0);
1932  StdRegions::ConstFactorMap factors = mkey.GetConstFactors();
1933  int edge = (int)factors[StdRegions::eFactorGaussEdge];
1934 
1935  coords[0] = (edge == 0 || edge == 3) ? -1.0 : 1.0;
1936 
1937  m_Ix = m_base[(edge + 1) % 2]->GetI(coords);
1938  returnval =
1939  MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0,m_Ix);
1940  }
1941  break;
1942  case StdRegions::ePreconLinearSpace:
1943  {
1944  NekDouble one = 1.0;
1945  MatrixKey helmkey(
1946  StdRegions::eHelmholtz, mkey.GetShapeType(), *this,
1947  mkey.GetConstFactors(), mkey.GetVarCoeffs());
1948  DNekScalBlkMatSharedPtr helmStatCond =
1949  GetLocStaticCondMatrix(helmkey);
1950  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1951  DNekMatSharedPtr R=BuildVertexMatrix(A);
1952 
1953  returnval =
1954  MemoryManager<DNekScalMat>::AllocateSharedPtr(one, R);
1955  }
1956  break;
1957  default:
1958  {
1959  NekDouble one = 1.0;
1960  DNekMatSharedPtr mat = GenMatrix(mkey);
1961 
1962  returnval =
1963  MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
1964  }
1965  break;
1966  }
1967 
1968  return returnval;
1969  }
1970 
1971 
1972  DNekScalBlkMatSharedPtr QuadExp::CreateStaticCondMatrix(
1973  const MatrixKey &mkey)
1974  {
1975  DNekScalBlkMatSharedPtr returnval;
1976 
1977  ASSERTL2(m_metricinfo->GetGtype()
1978  != SpatialDomains::eNoGeomType,
1979  "Geometric information is not set up");
1980 
1981  // set up block matrix system
1982  unsigned int nbdry = NumBndryCoeffs();
1983  unsigned int nint = (unsigned int)(m_ncoeffs - nbdry);
1984  unsigned int exp_size[] = {nbdry,nint};
1985  unsigned int nblks = 2;
1986  returnval = MemoryManager<DNekScalBlkMat>::
1987  AllocateSharedPtr(nblks,nblks,exp_size,exp_size);
1988  //Really need a constructor which takes Arrays
1989  NekDouble factor = 1.0;
1990 
1991  switch (mkey.GetMatrixType())
1992  {
1993  // this can only use stdregions statically condensed system
1994  // for mass matrix
1995  case StdRegions::eMass:
1996  if ((m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1997  ||(mkey.GetNVarCoeff()))
1998  {
1999  factor = 1.0;
2000  goto UseLocRegionsMatrix;
2001  }
2002  else
2003  {
2004  factor = (m_metricinfo->GetJac(GetPointsKeys()))[0];
2005  goto UseStdRegionsMatrix;
2006  }
2007  break;
2008  default: // use Deformed case for both
2009  // regular and deformed geometries
2010  factor = 1.0;
2011  goto UseLocRegionsMatrix;
2012  break;
2013  UseStdRegionsMatrix:
2014  {
2015  NekDouble invfactor = 1.0/factor;
2016  NekDouble one = 1.0;
2017  DNekBlkMatSharedPtr mat = GetStdStaticCondMatrix(mkey);
2018  DNekScalMatSharedPtr Atmp;
2019  DNekMatSharedPtr Asubmat;
2020 
2021  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::
2022  AllocateSharedPtr(factor,Asubmat = mat->GetBlock(0,0)));
2023  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::
2024  AllocateSharedPtr(one,Asubmat = mat->GetBlock(0,1)));
2025  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::
2026  AllocateSharedPtr(factor,Asubmat = mat->GetBlock(1,0)));
2027  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::
2028  AllocateSharedPtr(invfactor,Asubmat = mat->GetBlock(1,1)));
2029  }
2030  break;
2031  UseLocRegionsMatrix:
2032  {
2033  int i,j;
2034  NekDouble invfactor = 1.0/factor;
2035  NekDouble one = 1.0;
2036  DNekScalMat &mat = *GetLocMatrix(mkey);
2037  DNekMatSharedPtr A = MemoryManager<DNekMat>::
2038  AllocateSharedPtr(nbdry,nbdry);
2039  DNekMatSharedPtr B = MemoryManager<DNekMat>::
2040  AllocateSharedPtr(nbdry,nint);
2041  DNekMatSharedPtr C = MemoryManager<DNekMat>::
2042  AllocateSharedPtr(nint,nbdry);
2043  DNekMatSharedPtr D = MemoryManager<DNekMat>::
2044  AllocateSharedPtr(nint,nint);
2045 
2046  Array<OneD,unsigned int> bmap(nbdry);
2047  Array<OneD,unsigned int> imap(nint);
2048  GetBoundaryMap(bmap);
2049  GetInteriorMap(imap);
2050 
2051  for (i = 0; i < nbdry; ++i)
2052  {
2053  for(j = 0; j < nbdry; ++j)
2054  {
2055  (*A)(i,j) = mat(bmap[i],bmap[j]);
2056  }
2057 
2058  for(j = 0; j < nint; ++j)
2059  {
2060  (*B)(i,j) = mat(bmap[i],imap[j]);
2061  }
2062  }
2063 
2064  for (i = 0; i < nint; ++i)
2065  {
2066  for(j = 0; j < nbdry; ++j)
2067  {
2068  (*C)(i,j) = mat(imap[i],bmap[j]);
2069  }
2070 
2071  for(j = 0; j < nint; ++j)
2072  {
2073  (*D)(i,j) = mat(imap[i],imap[j]);
2074  }
2075  }
2076 
2077  // Calculate static condensed system
2078  if(nint)
2079  {
2080  D->Invert();
2081  (*B) = (*B)*(*D);
2082  (*A) = (*A) - (*B)*(*C);
2083  }
2084 
2085  DNekScalMatSharedPtr Atmp;
2086 
2087  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::
2088  AllocateSharedPtr(factor, A));
2089  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::
2090  AllocateSharedPtr(one, B));
2091  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::
2092  AllocateSharedPtr(factor, C));
2093  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::
2094  AllocateSharedPtr(invfactor, D));
2095 
2096  }
2097  }
2098  return returnval;
2099  }
2100 
2101 
2102  DNekScalMatSharedPtr QuadExp::v_GetLocMatrix(const MatrixKey &mkey)
2103  {
2104  return m_matrixManager[mkey];
2105  }
2106 
2107 
2108  DNekScalBlkMatSharedPtr QuadExp::v_GetLocStaticCondMatrix(
2109  const MatrixKey &mkey)
2110  {
2111  return m_staticCondMatrixManager[mkey];
2112  }
2113 
2114  void QuadExp::v_DropLocStaticCondMatrix(const MatrixKey &mkey)
2115  {
2116  m_staticCondMatrixManager.DeleteObject(mkey);
2117  }
2118 
2119 
2120  void QuadExp::v_MassMatrixOp(
2121  const Array<OneD, const NekDouble> &inarray,
2122  Array<OneD,NekDouble> &outarray,
2123  const StdRegions::StdMatrixKey &mkey)
2124  {
2125  StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
2126  }
2127 
2128 
2129  void QuadExp::v_LaplacianMatrixOp(
2130  const Array<OneD, const NekDouble> &inarray,
2131  Array<OneD,NekDouble> &outarray,
2132  const StdRegions::StdMatrixKey &mkey)
2133  {
2134  QuadExp::LaplacianMatrixOp_MatFree(inarray, outarray, mkey);
2135  }
2136 
2137 
2138  void QuadExp::v_LaplacianMatrixOp(
2139  const int k1,
2140  const int k2,
2141  const Array<OneD, const NekDouble> &inarray,
2142  Array<OneD,NekDouble> &outarray,
2143  const StdRegions::StdMatrixKey &mkey)
2144  {
2145  StdExpansion::LaplacianMatrixOp_MatFree(
2146  k1, k2, inarray, outarray, mkey);
2147  }
2148 
2149 
2150  void QuadExp::v_WeakDerivMatrixOp(
2151  const int i,
2152  const Array<OneD, const NekDouble> &inarray,
2153  Array<OneD,NekDouble> &outarray,
2154  const StdRegions::StdMatrixKey &mkey)
2155  {
2156  StdExpansion::WeakDerivMatrixOp_MatFree(i, inarray, outarray, mkey);
2157  }
2158 
2159 
2160  void QuadExp::v_WeakDirectionalDerivMatrixOp(
2161  const Array<OneD, const NekDouble> &inarray,
2162  Array<OneD,NekDouble> &outarray,
2163  const StdRegions::StdMatrixKey &mkey)
2164  {
2165  StdExpansion::WeakDirectionalDerivMatrixOp_MatFree(
2166  inarray, outarray, mkey);
2167  }
2168 
2169 
2170  void QuadExp::v_MassLevelCurvatureMatrixOp(
2171  const Array<OneD, const NekDouble> &inarray,
2172  Array<OneD,NekDouble> &outarray,
2173  const StdRegions::StdMatrixKey &mkey)
2174  {
2175  StdExpansion::MassLevelCurvatureMatrixOp_MatFree(
2176  inarray, outarray, mkey);
2177  }
2178 
2179 
2180  void QuadExp::v_HelmholtzMatrixOp(
2181  const Array<OneD, const NekDouble> &inarray,
2182  Array<OneD,NekDouble> &outarray,
2183  const StdRegions::StdMatrixKey &mkey)
2184  {
2185  QuadExp::HelmholtzMatrixOp_MatFree(inarray, outarray, mkey);
2186  }
2187 
2188 
2189  void QuadExp::v_GeneralMatrixOp_MatOp(
2190  const Array<OneD, const NekDouble> &inarray,
2191  Array<OneD,NekDouble> &outarray,
2192  const StdRegions::StdMatrixKey &mkey)
2193  {
2194  MatrixKey newkey(mkey);
2195  DNekScalMatSharedPtr mat = GetLocMatrix(newkey);
2196 
2197  if (inarray.get() == outarray.get())
2198  {
2199  Array<OneD,NekDouble> tmp(m_ncoeffs);
2200  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
2201 
2202  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs, mat->Scale(),
2203  (mat->GetOwnedMatrix())->GetPtr().get(), m_ncoeffs,
2204  tmp.get(), 1, 0.0, outarray.get(), 1);
2205  }
2206  else
2207  {
2208  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),
2209  (mat->GetOwnedMatrix())->GetPtr().get(), m_ncoeffs,
2210  inarray.get(), 1, 0.0, outarray.get(), 1);
2211  }
2212  }
2213 
2214  void QuadExp::v_ReduceOrderCoeffs(
2215  int numMin,
2216  const Array<OneD, const NekDouble> &inarray,
2217  Array<OneD, NekDouble> &outarray)
2218  {
2219  int n_coeffs = inarray.num_elements();
2220 
2221  Array<OneD, NekDouble> coeff (n_coeffs);
2222  Array<OneD, NekDouble> coeff_tmp(n_coeffs, 0.0);
2223  Array<OneD, NekDouble> tmp, tmp2;
2224 
2225  int nmodes0 = m_base[0]->GetNumModes();
2226  int nmodes1 = m_base[1]->GetNumModes();
2227  int numMax = nmodes0;
2228 
2229  Vmath::Vcopy(n_coeffs,inarray,1,coeff_tmp,1);
2230 
2231  const LibUtilities::PointsKey Pkey0(
2232  nmodes0, LibUtilities::eGaussLobattoLegendre);
2233  const LibUtilities::PointsKey Pkey1(
2234  nmodes1, LibUtilities::eGaussLobattoLegendre);
2235  LibUtilities::BasisKey b0(
2236  m_base[0]->GetBasisType(), nmodes0, Pkey0);
2237  LibUtilities::BasisKey b1(
2238  m_base[1]->GetBasisType(), nmodes1, Pkey1);
2239  LibUtilities::BasisKey bortho0(
2240  LibUtilities::eOrtho_A, nmodes0, Pkey0);
2241  LibUtilities::BasisKey bortho1(
2242  LibUtilities::eOrtho_A, nmodes1, Pkey1);
2243 
2244  LibUtilities::InterpCoeff2D(
2245  b0, b1, coeff_tmp, bortho0, bortho1, coeff);
2246 
2247  Vmath::Zero(n_coeffs, coeff_tmp, 1);
2248 
2249  int cnt = 0;
2250  for (int i = 0; i < numMin+1; ++i)
2251  {
2252  Vmath::Vcopy(numMin,
2253  tmp = coeff+cnt,1,
2254  tmp2 = coeff_tmp+cnt,1);
2255 
2256  cnt = i*numMax;
2257  }
2258 
2259  LibUtilities::InterpCoeff2D(
2260  bortho0, bortho1, coeff_tmp,
2261  b0, b1, outarray);
2262  }
2263 
2264  void QuadExp::v_LaplacianMatrixOp_MatFree_Kernel(
2265  const Array<OneD, const NekDouble> &inarray,
2266  Array<OneD, NekDouble> &outarray,
2267  Array<OneD, NekDouble> &wsp)
2268  {
2269  if (m_metrics.count(eMetricLaplacian00) == 0)
2270  {
2271  ComputeLaplacianMetric();
2272  }
2273 
2274  int nquad0 = m_base[0]->GetNumPoints();
2275  int nquad1 = m_base[1]->GetNumPoints();
2276  int nqtot = nquad0*nquad1;
2277  int nmodes0 = m_base[0]->GetNumModes();
2278  int nmodes1 = m_base[1]->GetNumModes();
2279  int wspsize = max(max(max(nqtot,m_ncoeffs),nquad1*nmodes0),nquad0*nmodes1);
2280 
2281  ASSERTL1(wsp.num_elements() >= 3*wspsize,
2282  "Workspace is of insufficient size.");
2283 
2284  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
2285  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
2286  const Array<OneD, const NekDouble>& dbase0 = m_base[0]->GetDbdata();
2287  const Array<OneD, const NekDouble>& dbase1 = m_base[1]->GetDbdata();
2288  const Array<OneD, const NekDouble>& metric00 = m_metrics[eMetricLaplacian00];
2289  const Array<OneD, const NekDouble>& metric01 = m_metrics[eMetricLaplacian01];
2290  const Array<OneD, const NekDouble>& metric11 = m_metrics[eMetricLaplacian11];
2291 
2292  // Allocate temporary storage
2293  Array<OneD,NekDouble> wsp0(wsp);
2294  Array<OneD,NekDouble> wsp1(wsp+wspsize);
2295  Array<OneD,NekDouble> wsp2(wsp+2*wspsize);
2296 
2297  StdExpansion2D::PhysTensorDeriv(inarray,wsp1,wsp2);
2298 
2299  // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
2300  // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
2301  // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
2302  // especially for this purpose
2303  Vmath::Vvtvvtp(nqtot,&metric00[0],1,&wsp1[0],1,&metric01[0],1,&wsp2[0],1,&wsp0[0],1);
2304  Vmath::Vvtvvtp(nqtot,&metric01[0],1,&wsp1[0],1,&metric11[0],1,&wsp2[0],1,&wsp2[0],1);
2305 
2306  // outarray = m = (D_xi1 * B)^T * k
2307  // wsp1 = n = (D_xi2 * B)^T * l
2308  IProductWRTBase_SumFacKernel(dbase0,base1,wsp0,outarray,wsp1,false,true);
2309  IProductWRTBase_SumFacKernel(base0,dbase1,wsp2,wsp1, wsp0,true,false);
2310 
2311  // outarray = outarray + wsp1
2312  // = L * u_hat
2313  Vmath::Vadd(m_ncoeffs,wsp1.get(),1,outarray.get(),1,outarray.get(),1);
2314  }
2315 
2316  void QuadExp::v_ComputeLaplacianMetric()
2317  {
2318  if (m_metrics.count(eMetricQuadrature) == 0)
2319  {
2320  ComputeQuadratureMetric();
2321  }
2322 
2323  const SpatialDomains::GeomType type = m_metricinfo->GetGtype();
2324  const unsigned int nqtot = GetTotPoints();
2325  const unsigned int dim = 2;
2326  const MetricType m[3][3] = { {eMetricLaplacian00, eMetricLaplacian01, eMetricLaplacian02},
2327  {eMetricLaplacian01, eMetricLaplacian11, eMetricLaplacian12},
2328  {eMetricLaplacian02, eMetricLaplacian12, eMetricLaplacian22}
2329  };
2330 
2331  const Array<TwoD, const NekDouble> gmat =
2332  m_metricinfo->GetGmat(GetPointsKeys());
2333  for (unsigned int i = 0; i < dim; ++i)
2334  {
2335  for (unsigned int j = i; j < dim; ++j)
2336  {
2337  m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
2338  if (type == SpatialDomains::eDeformed)
2339  {
2340  Vmath::Vcopy(nqtot, &gmat[i*dim+j][0], 1,
2341  &m_metrics[m[i][j]][0], 1);
2342  }
2343  else
2344  {
2345  Vmath::Fill(nqtot, gmat[i*dim+j][0],
2346  &m_metrics[m[i][j]][0], 1);
2347  }
2348  MultiplyByQuadratureMetric(m_metrics[m[i][j]],
2349  m_metrics[m[i][j]]);
2350 
2351  }
2352  }
2353  }
2354 
2355  void QuadExp::v_SVVLaplacianFilter(
2356  Array<OneD, NekDouble> &array,
2357  const StdRegions::StdMatrixKey &mkey)
2358  {
2359  int nq = GetTotPoints();
2360 
2361  // Calculate sqrt of the Jacobian
2362  Array<OneD, const NekDouble> jac =
2363  m_metricinfo->GetJac(GetPointsKeys());
2364  Array<OneD, NekDouble> sqrt_jac(nq);
2365  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
2366  {
2367  Vmath::Vsqrt(nq,jac,1,sqrt_jac,1);
2368  }
2369  else
2370  {
2371  Vmath::Fill(nq,sqrt(jac[0]),sqrt_jac,1);
2372  }
2373 
2374  // Multiply array by sqrt(Jac)
2375  Vmath::Vmul(nq,sqrt_jac,1,array,1,array,1);
2376 
2377  // Apply std region filter
2378  StdQuadExp::v_SVVLaplacianFilter( array, mkey);
2379 
2380  // Divide by sqrt(Jac)
2381  Vmath::Vdiv(nq,array,1,sqrt_jac,1,array,1);
2382  }
2383 
2384  }//end of namespace
2385 }//end of namespace
Nektar::StdRegions::StdExpansion::GetPointsKeys
const LibUtilities::PointsKeyVector GetPointsKeys() const
Definition: StdExpansion.h:1323
Nektar::LocalRegions::QuadExp::v_ReduceOrderCoeffs
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: QuadExp.cpp:2214
Nektar::StdRegions::StdExpansion::DetShapeType
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:470
Nektar::StdRegions::StdMatrixKey::GetConstFactor
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:122
Nektar::StdRegions::StdExpansion::GenMatrix
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:1054
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Nektar::StdRegions::StdExpansion::GetCartesianEorient
StdRegions::Orientation GetCartesianEorient(int edge)
Definition: StdExpansion.h:778
Nektar::LocalRegions::QuadExp::v_GetCartesianEorient
virtual StdRegions::Orientation v_GetCartesianEorient(int edge)
Definition: QuadExp.cpp:1608
Nektar::StdRegions::StdMatrixKey::GetConstFactors
const ConstFactorMap & GetConstFactors() const
Definition: StdMatrixKey.h:142
Nektar::LocalRegions::QuadExp::v_MassMatrixOp
virtual void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: QuadExp.cpp:2120
Nektar::LocalRegions::QuadExp::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: QuadExp.cpp:576
Nektar::StdRegions::StdMatrixKey::GetVarCoeffs
const VarCoeffMap & GetVarCoeffs() const
Definition: StdMatrixKey.h:168
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:242
Nektar::StdRegions::StdMatrixKey::GetMatrixType
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:82
Nektar::LocalRegions::QuadExp::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*base1 and put into outarray...
Definition: QuadExp.cpp:409
Nektar::StdRegions::StdExpansion::MassMatrixOp
void MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:976
Nektar::LocalRegions::QuadExp::v_GetBasis
virtual const LibUtilities::BasisSharedPtr & v_GetBasis(int dir) const
Definition: QuadExp.cpp:1614
Nektar::LocalRegions::QuadExp::v_FwdTrans_BndConstrained
virtual void v_FwdTrans_BndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: QuadExp.cpp:286
Nektar::NullNekDouble1DArray
static Array< OneD, NekDouble > NullNekDouble1DArray
Definition: SharedArray.hpp:785
Nektar::LocalRegions::QuadExp::m_staticCondMatrixManager
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
Definition: QuadExp.h:290
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:27
Nektar::StdRegions::StdExpansion::IProductWRTDerivBase_SumFac
void IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:1463
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:408
Nektar::StdRegions::StdExpansion::MultiplyByQuadratureMetric
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:947
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:102
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:635
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
Nektar::LocalRegions::QuadExp::v_DropLocStaticCondMatrix
void v_DropLocStaticCondMatrix(const MatrixKey &mkey)
Definition: QuadExp.cpp:2114
Nektar::SpatialDomains::QuadGeomSharedPtr
boost::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: QuadGeom.h:54
Nektar::LocalRegions::QuadExp::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: QuadExp.cpp:260
VmathArray.hpp
Nektar::StdRegions::StdExpansion::LaplacianMatrixOp_MatFree
void LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:1485
Nektar::LocalRegions::QuadExp::v_GetNumPoints
virtual int v_GetNumPoints(const int dir) const
Definition: QuadExp.cpp:1621
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:229
Nektar::LocalRegions::QuadExp::v_GetMetricInfo
virtual const SpatialDomains::GeomFactorsSharedPtr & v_GetMetricInfo() const
Definition: QuadExp.cpp:1491
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:485
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:442
Nektar::LibUtilities::eModified_A
Principle Modified Functions .
Definition: BasisType.h:49
Nektar::StdRegions::StdExpansion::NormVectorIProductWRTBase
void NormVectorIProductWRTBase(const Array< OneD, const NekDouble > &Fx, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:733
Nektar::LocalRegions::Expansion2D::m_edgeNormals
std::map< int, StdRegions::NormalVector > m_edgeNormals
Definition: Expansion2D.h:135
Nektar::LocalRegions::Expansion::m_metricinfo
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
Definition: Expansion.h:131
Nektar::LibUtilities::eGauss_Lagrange
Lagrange Polynomials using the Gauss points .
Definition: BasisType.h:54
Nektar::StdRegions::StdExpansion::GetEorient
StdRegions::Orientation GetEorient(int edge)
Definition: StdExpansion.h:768
Nektar::LocalRegions::QuadExp::CreateStaticCondMatrix
DNekScalBlkMatSharedPtr CreateStaticCondMatrix(const MatrixKey &mkey)
Definition: QuadExp.cpp:1972
std
STL namespace.
Nektar::LocalRegions::QuadExp::v_GetEdgeInterpVals
virtual void v_GetEdgeInterpVals(const int edge, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: QuadExp.cpp:827
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:271
Nektar::LocalRegions::Expansion::BuildVertexMatrix
DNekMatSharedPtr BuildVertexMatrix(const DNekScalMatSharedPtr &r_bnd)
Definition: Expansion.cpp:98
Nektar::StdRegions::StdMatrixKey::GetShapeType
LibUtilities::ShapeType GetShapeType() const
Definition: StdMatrixKey.h:87
Nektar::StdRegions::ConstFactorMap
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:252
Nektar::LocalRegions::QuadExp::v_GetEdgePhysMap
virtual void v_GetEdgePhysMap(const int edge, Array< OneD, int > &outarray)
Definition: QuadExp.cpp:893
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:241
Nektar::LocalRegions::Expansion::m_geom
SpatialDomains::GeometrySharedPtr m_geom
Definition: Expansion.h:130
Nektar::StdRegions::StdExpansion::GetLocStaticCondMatrix
DNekScalBlkMatSharedPtr GetLocStaticCondMatrix(const LocalRegions::MatrixKey &mkey)
Definition: StdExpansion.h:753
Nektar::DNekMatSharedPtr
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
Nektar::LocalRegions::QuadExp::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: QuadExp.cpp:203
Nektar::StdRegions::StdExpansion::GetStdMatrix
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:706
Nektar::DNekScalMatSharedPtr
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Definition: NekMatrixFwd.hpp:73
Nektar::LibUtilities::eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:49
Nektar::LocalRegions::QuadExp::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: QuadExp.cpp:752
Nektar::StdRegions::StdExpansion::IProductWRTBase_SumFac
void IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Definition: StdExpansion.h:1410
Nektar::LocalRegions::QuadExp::v_GetLocMatrix
virtual DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey)
Definition: QuadExp.cpp:2102
Nektar::LocalRegions::QuadExp::v_GeneralMatrixOp_MatOp
virtual void v_GeneralMatrixOp_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: QuadExp.cpp:2189
Nektar::StdRegions::StdMatrixKey::ConstFactorExists
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:131
Nektar::StdRegions::StdExpansion::GetTotPoints
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:141
Nektar::LocalRegions::QuadExp::v_HelmholtzMatrixOp
virtual void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: QuadExp.cpp:2180
Nektar::LibUtilities::Interp2D
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis...
Definition: Interp.cpp:116
Nektar::LocalRegions::QuadExp::v_GetEorient
virtual StdRegions::Orientation v_GetEorient(int edge)
Definition: QuadExp.cpp:1602
Nektar::LocalRegions::QuadExp::v_GetLinStdExp
virtual StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const
Definition: QuadExp.cpp:620
Nektar::LocalRegions::QuadExp::v_GetCoord
virtual void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords)
Definition: QuadExp.cpp:640
Nektar::StdRegions::StdExpansion::GetStdStaticCondMatrix
DNekBlkMatSharedPtr GetStdStaticCondMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:711
Nektar::LocalRegions::SegExpSharedPtr
boost::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:270
Nektar::LocalRegions::Expansion2D::GetGeom2D
SpatialDomains::Geometry2DSharedPtr GetGeom2D() const
Definition: Expansion2D.h:269
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1085
Nektar::LocalRegions::QuadExp::v_GetCoords
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
Definition: QuadExp.cpp:631
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:213
Nektar::LocalRegions::QuadExp::v_Integral
virtual NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
Integrates the specified function over the domain.
Definition: QuadExp.cpp:84
Nektar::LocalRegions::QuadExp::v_IProductWRTDerivBase_MatOp
virtual void v_IProductWRTDerivBase_MatOp(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: QuadExp.cpp:535
Nektar::LocalRegions::QuadExp::m_matrixManager
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
Definition: QuadExp.h:289
Nektar::LocalRegions::QuadExp::v_IProductWRTBase_SumFac
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Definition: QuadExp.cpp:433
Nektar::LocalRegions::QuadExp::v_GetLocStaticCondMatrix
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey)
Definition: QuadExp.cpp:2108
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:224
Nektar::LocalRegions::QuadExp::v_IProductWRTDerivBase
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: QuadExp.cpp:424
Nektar::DNekScalBlkMatSharedPtr
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
Nektar::StdRegions::StdExpansion::GetInteriorMap
void GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:819
Nektar::LocalRegions::QuadExp::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: QuadExp.cpp:110
Nektar::LibUtilities::eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:46
Nektar::LocalRegions::QuadExp::v_IProductWRTBase_MatOp
virtual void v_IProductWRTBase_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: QuadExp.cpp:463
Nektar::LocalRegions::QuadExp::CreateMatrix
DNekScalMatSharedPtr CreateMatrix(const MatrixKey &mkey)
Definition: QuadExp.cpp:1658
Nektar::Transpose
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Definition: BlockMatrix.hpp:219
Nektar::LocalRegions::QuadExp::v_LaplacianMatrixOp
virtual void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: QuadExp.cpp:2129
Expansion3D.h
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:58
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::LocalRegions::Expansion2D::GetLeftAdjacentElementExp
Expansion3DSharedPtr GetLeftAdjacentElementExp() const
Definition: Expansion2D.h:223
SegExp.h
Interp.h
QuadExp.h
Nektar::LocalRegions::QuadExp::v_MassLevelCurvatureMatrixOp
virtual void v_MassLevelCurvatureMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: QuadExp.cpp:2170
Nektar::LocalRegions::Expansion2D::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: Expansion2D.cpp:641
Vmath::Ddot
T Ddot(int n, const Array< OneD, const T > &w, const int incw, const Array< OneD, const T > &x, const int incx, const Array< OneD, const int > &y, const int incy)
Definition: VmathArray.hpp:436
Nektar::LocalRegions::QuadExp::v_WeakDerivMatrixOp
virtual void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: QuadExp.cpp:2150
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:165
Nektar::LocalRegions::QuadExp::v_GetStdExp
virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const
Definition: QuadExp.cpp:612
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:343
Nektar::DNekBlkMatSharedPtr
boost::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
Definition: NekTypeDefs.hpp:72
Nektar::LocalRegions::Expansion::GetLocMatrix
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition: Expansion.cpp:85
Nektar::LocalRegions::QuadExp::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: QuadExp.cpp:2264
Vmath.hpp
Nektar::LocalRegions::QuadExp::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: QuadExp.cpp:1626
Nektar::LocalRegions::Expansion::GetGeom
SpatialDomains::GeometrySharedPtr GetGeom() const
Definition: Expansion.cpp:161
Nektar::SpatialDomains::GeomFactorsSharedPtr
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
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:851
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:537
Nektar::LocalRegions::QuadExp::v_IProductWRTDerivBase_SumFac
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: QuadExp.cpp:479
ASSERTL2
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:250
Nektar::SpatialDomains::eRegular
Geometry is straight-sided with constant geometric factors.
Definition: SpatialDomains.hpp:55
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:54
Nektar::StdRegions::StdExpansion::HelmholtzMatrixOp_MatFree
void HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:1527
Nektar::LibUtilities::PointsKey::GetNumPoints
unsigned int GetNumPoints() const
Definition: Points.h:106
Nektar::LocalRegions::QuadExp::v_GetEdgeQFactors
virtual void v_GetEdgeQFactors(const int edge, Array< OneD, NekDouble > &outarray)
Definition: QuadExp.cpp:941
Nektar::StdRegions::StdQuadExpSharedPtr
boost::shared_ptr< StdQuadExp > StdQuadExpSharedPtr
Definition: StdQuadExp.h:262
Nektar::StdRegions::StdExpansion::BwdTrans
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Backward transformation from coefficient space to physical space...
Definition: StdExpansion.h:531
Nektar::LocalRegions::QuadExp::v_WeakDirectionalDerivMatrixOp
virtual void v_WeakDirectionalDerivMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: QuadExp.cpp:2160
Nektar::LocalRegions::QuadExp::v_ComputeLaplacianMetric
virtual void v_ComputeLaplacianMetric()
Definition: QuadExp.cpp:2316
Nektar::StdRegions::StdExpansion::GetNcoeffs
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:131
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:591
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:200
InterpCoeff.h
Nektar::SpatialDomains::GeomType
GeomType
Indicates the type of element geometry.
Definition: SpatialDomains.hpp:52
Nektar::LocalRegions::QuadExp::v_StdPhysEvaluate
virtual NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals)
Definition: QuadExp.cpp:663
Nektar::LibUtilities::BasisSharedPtr
boost::shared_ptr< Basis > BasisSharedPtr
Definition: FoundationsFwd.hpp:80
Nektar::LibUtilities::eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:53
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
Nektar::StdRegions::StdExpansionSharedPtr
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1878
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:216
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition: StdExpansion.h:1418
Nektar::LibUtilities::InterpCoeff2D
void InterpCoeff2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
Definition: InterpCoeff.cpp:73
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
Nektar::SpatialDomains::eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:57
Nektar::LocalRegions::QuadExp::v_ComputeEdgeNormal
virtual void v_ComputeEdgeNormal(const int edge)
Definition: QuadExp.cpp:1201
Nektar::LocalRegions::QuadExp::v_CreateStdMatrix
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: QuadExp.cpp:1647
Nektar::StdRegions::StdExpansion::GetBoundaryMap
void GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:814
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
Nektar::LibUtilities::eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:52
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:299
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:183
Nektar::LocalRegions::QuadExp::v_SVVLaplacianFilter
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey)
Definition: QuadExp.cpp:2355
Nektar::LocalRegions::QuadExp::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: QuadExp.cpp:687
Nektar::LocalRegions::QuadExp::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: QuadExp.cpp:671
Nektar::LocalRegions::QuadExp::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: QuadExp.cpp:1503
Nektar::StdRegions::StdExpansion::FwdTrans
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space...
Definition: StdExpansion.h:1910
Generated on Sat Nov 25 2017 16:18:16 for Nektar++ by   doxygen 1.8.8