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  // Assume that input is also Gll_Lagrange but no way to check;
1567  LibUtilities::PointsKey
1568  p0(nummodes[0], LibUtilities::eGaussLobattoLegendre);
1569  LibUtilities::PointsKey
1570  p1(nummodes[1], LibUtilities::eGaussLobattoLegendre);
1571  LibUtilities::Interp2D(p0, p1, data,
1572  m_base[0]->GetPointsKey(),
1573  m_base[1]->GetPointsKey(),
1574  coeffs);
1575  }
1576  break;
1577  case LibUtilities::eGauss_Lagrange:
1578  {
1579  // Assume that input is also Gll_Lagrange but no way to check;
1580  LibUtilities::PointsKey
1581  p0(nummodes[0],LibUtilities::eGaussGaussLegendre);
1582  LibUtilities::PointsKey
1583  p1(nummodes[1],LibUtilities::eGaussGaussLegendre);
1584  LibUtilities::Interp2D(p0, p1, data,
1585  m_base[0]->GetPointsKey(),
1586  m_base[1]->GetPointsKey(),
1587  coeffs);
1588  }
1589  break;
1590  default:
1591  ASSERTL0(false,
1592  "basis is either not set up or not hierarchicial");
1593  }
1594  }
1595 
1596 
1597  StdRegions::Orientation QuadExp::v_GetEorient(int edge)
1598  {
1599  return m_geom->GetEorient(edge);
1600  }
1601 
1602 
1603  StdRegions::Orientation QuadExp::v_GetCartesianEorient(int edge)
1604  {
1605  return GetGeom2D()->GetCartesianEorient(edge);
1606  }
1607 
1608 
1609  const LibUtilities::BasisSharedPtr& QuadExp::v_GetBasis(int dir) const
1610  {
1611  ASSERTL1(dir >= 0 &&dir <= 1, "input dir is out of range");
1612  return m_base[dir];
1613  }
1614 
1615 
1616  int QuadExp::v_GetNumPoints(const int dir) const
1617  {
1618  return GetNumPoints(dir);
1619  }
1620 
1621  DNekMatSharedPtr QuadExp::v_GenMatrix(
1622  const StdRegions::StdMatrixKey &mkey)
1623  {
1624  DNekMatSharedPtr returnval;
1625  switch (mkey.GetMatrixType())
1626  {
1627  case StdRegions::eHybridDGHelmholtz:
1628  case StdRegions::eHybridDGLamToU:
1629  case StdRegions::eHybridDGLamToQ0:
1630  case StdRegions::eHybridDGLamToQ1:
1631  case StdRegions::eHybridDGLamToQ2:
1632  case StdRegions::eHybridDGHelmBndLam:
1633  case StdRegions::eInvLaplacianWithUnityMean:
1634  returnval = Expansion2D::v_GenMatrix(mkey);
1635  break;
1636  default:
1637  returnval = StdQuadExp::v_GenMatrix(mkey);
1638  }
1639  return returnval;
1640  }
1641 
1642  DNekMatSharedPtr QuadExp::v_CreateStdMatrix(
1643  const StdRegions::StdMatrixKey &mkey)
1644  {
1645  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
1646  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
1647  StdRegions::StdQuadExpSharedPtr tmp =
1648  MemoryManager<StdQuadExp>::AllocateSharedPtr(bkey0,bkey1);
1649  return tmp->GetStdMatrix(mkey);
1650  }
1651 
1652 
1653  DNekScalMatSharedPtr QuadExp::CreateMatrix(const MatrixKey &mkey)
1654  {
1655  DNekScalMatSharedPtr returnval;
1656  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
1657 
1658  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,
1659  "Geometric information is not set up");
1660 
1661  switch (mkey.GetMatrixType())
1662  {
1663  case StdRegions::eMass:
1664  {
1665  if ((m_metricinfo->GetGtype() ==
1666  SpatialDomains::eDeformed) || (mkey.GetNVarCoeff()))
1667  {
1668  NekDouble one = 1.0;
1669  DNekMatSharedPtr mat = GenMatrix(mkey);
1670  returnval = MemoryManager<DNekScalMat>::
1671  AllocateSharedPtr(one,mat);
1672  }
1673  else
1674  {
1675  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1676  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1677  returnval = MemoryManager<DNekScalMat>::
1678  AllocateSharedPtr(jac,mat);
1679  }
1680  }
1681  break;
1682  case StdRegions::eInvMass:
1683  {
1684  if ((m_metricinfo->GetGtype() ==
1685  SpatialDomains::eDeformed) || (mkey.GetNVarCoeff()))
1686  {
1687  NekDouble one = 1.0;
1688  StdRegions::StdMatrixKey masskey(
1689  StdRegions::eMass, DetShapeType(), *this);
1690  DNekMatSharedPtr mat = GenMatrix(masskey);
1691  mat->Invert();
1692 
1693  returnval = MemoryManager<DNekScalMat>::
1694  AllocateSharedPtr(one,mat);
1695  }
1696  else
1697  {
1698  NekDouble fac = 1.0/(m_metricinfo->GetJac(ptsKeys))[0];
1699  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1700  returnval = MemoryManager<DNekScalMat>::
1701  AllocateSharedPtr(fac,mat);
1702  }
1703  }
1704  break;
1705  case StdRegions::eWeakDeriv0:
1706  case StdRegions::eWeakDeriv1:
1707  case StdRegions::eWeakDeriv2:
1708  {
1709  if((m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1710  || (mkey.GetNVarCoeff()))
1711  {
1712  NekDouble one = 1.0;
1713  DNekMatSharedPtr mat = GenMatrix(mkey);
1714 
1715  returnval = MemoryManager<DNekScalMat>::
1716  AllocateSharedPtr(one,mat);
1717  }
1718  else
1719  {
1720  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1721  Array<TwoD, const NekDouble> df =
1722  m_metricinfo->GetDerivFactors(ptsKeys);
1723  int dir = 0;
1724 
1725  switch(mkey.GetMatrixType())
1726  {
1727  case StdRegions::eWeakDeriv0:
1728  dir = 0;
1729  break;
1730  case StdRegions::eWeakDeriv1:
1731  dir = 1;
1732  break;
1733  case StdRegions::eWeakDeriv2:
1734  dir = 2;
1735  break;
1736  default:
1737  break;
1738  }
1739 
1740  MatrixKey deriv0key(StdRegions::eWeakDeriv0,
1741  mkey.GetShapeType(), *this);
1742  MatrixKey deriv1key(StdRegions::eWeakDeriv1,
1743  mkey.GetShapeType(), *this);
1744 
1745  DNekMat &deriv0 = *GetStdMatrix(deriv0key);
1746  DNekMat &deriv1 = *GetStdMatrix(deriv1key);
1747 
1748  int rows = deriv0.GetRows();
1749  int cols = deriv1.GetColumns();
1750 
1751  DNekMatSharedPtr WeakDeriv = MemoryManager<DNekMat>::
1752  AllocateSharedPtr(rows,cols);
1753  (*WeakDeriv) = df[2*dir][0]*deriv0 +
1754  df[2*dir+1][0]*deriv1;
1755  returnval = MemoryManager<DNekScalMat>::
1756  AllocateSharedPtr(jac,WeakDeriv);
1757  }
1758  }
1759  break;
1760  case StdRegions::eLaplacian:
1761  {
1762  if( (m_metricinfo->GetGtype() ==
1763  SpatialDomains::eDeformed) || (mkey.GetNVarCoeff() > 0)
1764  || (mkey.ConstFactorExists
1765  (StdRegions::eFactorSVVCutoffRatio)))
1766  {
1767  NekDouble one = 1.0;
1768  DNekMatSharedPtr mat = GenMatrix(mkey);
1769 
1770  returnval = MemoryManager<DNekScalMat>::
1771  AllocateSharedPtr(one,mat);
1772  }
1773  else
1774  {
1775  MatrixKey lap00key(StdRegions::eLaplacian00,
1776  mkey.GetShapeType(), *this);
1777  MatrixKey lap01key(StdRegions::eLaplacian01,
1778  mkey.GetShapeType(), *this);
1779  MatrixKey lap11key(StdRegions::eLaplacian11,
1780  mkey.GetShapeType(), *this);
1781 
1782  DNekMat &lap00 = *GetStdMatrix(lap00key);
1783  DNekMat &lap01 = *GetStdMatrix(lap01key);
1784  DNekMat &lap11 = *GetStdMatrix(lap11key);
1785 
1786  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1787  Array<TwoD, const NekDouble>
1788  gmat = m_metricinfo->GetGmat(ptsKeys);
1789 
1790  int rows = lap00.GetRows();
1791  int cols = lap00.GetColumns();
1792 
1793  DNekMatSharedPtr lap =
1794  MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols);
1795 
1796  (*lap) = gmat[0][0] * lap00 +
1797  gmat[1][0] * (lap01 + Transpose(lap01)) +
1798  gmat[3][0] * lap11;
1799 
1800  returnval = MemoryManager<DNekScalMat>::
1801  AllocateSharedPtr(jac,lap);
1802  }
1803  }
1804  break;
1805  case StdRegions::eInvLaplacianWithUnityMean:
1806  {
1807  DNekMatSharedPtr mat = GenMatrix(mkey);
1808  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0,mat);
1809  }
1810  break;
1811  case StdRegions::eHelmholtz:
1812  {
1813  NekDouble lambda =
1814  mkey.GetConstFactor(StdRegions::eFactorLambda);
1815 
1816  MatrixKey masskey(mkey, StdRegions::eMass);
1817  DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
1818 
1819  MatrixKey lapkey(mkey, StdRegions::eLaplacian);
1820  DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
1821 
1822  int rows = LapMat.GetRows();
1823  int cols = LapMat.GetColumns();
1824 
1825  DNekMatSharedPtr helm = MemoryManager<DNekMat>::
1826  AllocateSharedPtr(rows,cols);
1827 
1828  NekDouble one = 1.0;
1829  (*helm) = LapMat + lambda*MassMat;
1830 
1831  returnval =
1832  MemoryManager<DNekScalMat>::AllocateSharedPtr(one,helm);
1833  }
1834  break;
1835  case StdRegions::eIProductWRTBase:
1836  {
1837  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1838  {
1839  NekDouble one = 1.0;
1840  DNekMatSharedPtr mat = GenMatrix(mkey);
1841  returnval = MemoryManager<DNekScalMat>::
1842  AllocateSharedPtr(one,mat);
1843  }
1844  else
1845  {
1846  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1847  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1848  returnval = MemoryManager<DNekScalMat>::
1849  AllocateSharedPtr(jac,mat);
1850  }
1851  }
1852  break;
1853  case StdRegions::eIProductWRTDerivBase0:
1854  case StdRegions::eIProductWRTDerivBase1:
1855  case StdRegions::eIProductWRTDerivBase2:
1856  {
1857  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1858  {
1859  NekDouble one = 1.0;
1860  DNekMatSharedPtr mat = GenMatrix(mkey);
1861  returnval = MemoryManager<DNekScalMat>::
1862  AllocateSharedPtr(one,mat);
1863  }
1864  else
1865  {
1866  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1867  const Array<TwoD, const NekDouble>& df =
1868  m_metricinfo->GetDerivFactors(ptsKeys);
1869  int dir = 0;
1870 
1871  switch(mkey.GetMatrixType())
1872  {
1873  case StdRegions::eIProductWRTDerivBase0:
1874  dir = 0;
1875  break;
1876  case StdRegions::eIProductWRTDerivBase1:
1877  dir = 1;
1878  break;
1879  case StdRegions::eIProductWRTDerivBase2:
1880  dir = 2;
1881  break;
1882  default:
1883  break;
1884  }
1885 
1886  MatrixKey iProdDeriv0Key(
1887  StdRegions::eIProductWRTDerivBase0,
1888  mkey.GetShapeType(), *this);
1889  MatrixKey iProdDeriv1Key(
1890  StdRegions::eIProductWRTDerivBase1,
1891  mkey.GetShapeType(), *this);
1892 
1893  DNekMat &stdiprod0 = *GetStdMatrix(iProdDeriv0Key);
1894  DNekMat &stdiprod1 = *GetStdMatrix(iProdDeriv0Key);
1895 
1896  int rows = stdiprod0.GetRows();
1897  int cols = stdiprod1.GetColumns();
1898 
1899  DNekMatSharedPtr mat = MemoryManager<DNekMat>::
1900  AllocateSharedPtr(rows,cols);
1901  (*mat) = df[2*dir][0]*stdiprod0 +
1902  df[2*dir+1][0]*stdiprod1;
1903 
1904  returnval = MemoryManager<DNekScalMat>::
1905  AllocateSharedPtr(jac,mat);
1906  }
1907  }
1908  break;
1909  case StdRegions::eInvHybridDGHelmholtz:
1910  {
1911  NekDouble one = 1.0;
1912 
1913  MatrixKey hkey(StdRegions::eHybridDGHelmholtz,
1914  DetShapeType(), *this,
1915  mkey.GetConstFactors(), mkey.GetVarCoeffs());
1916  DNekMatSharedPtr mat = GenMatrix(hkey);
1917 
1918  mat->Invert();
1919  returnval =
1920  MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1921  }
1922  break;
1923  case StdRegions::eInterpGauss:
1924  {
1925  DNekMatSharedPtr m_Ix;
1926  Array<OneD, NekDouble> coords(1, 0.0);
1927  StdRegions::ConstFactorMap factors = mkey.GetConstFactors();
1928  int edge = (int)factors[StdRegions::eFactorGaussEdge];
1929 
1930  coords[0] = (edge == 0 || edge == 3) ? -1.0 : 1.0;
1931 
1932  m_Ix = m_base[(edge + 1) % 2]->GetI(coords);
1933  returnval =
1934  MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0,m_Ix);
1935  }
1936  break;
1937  case StdRegions::ePreconLinearSpace:
1938  {
1939  NekDouble one = 1.0;
1940  MatrixKey helmkey(
1941  StdRegions::eHelmholtz, mkey.GetShapeType(), *this,
1942  mkey.GetConstFactors(), mkey.GetVarCoeffs());
1943  DNekScalBlkMatSharedPtr helmStatCond =
1944  GetLocStaticCondMatrix(helmkey);
1945  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1946  DNekMatSharedPtr R=BuildVertexMatrix(A);
1947 
1948  returnval =
1949  MemoryManager<DNekScalMat>::AllocateSharedPtr(one, R);
1950  }
1951  break;
1952  default:
1953  {
1954  NekDouble one = 1.0;
1955  DNekMatSharedPtr mat = GenMatrix(mkey);
1956 
1957  returnval =
1958  MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
1959  }
1960  break;
1961  }
1962 
1963  return returnval;
1964  }
1965 
1966 
1967  DNekScalBlkMatSharedPtr QuadExp::CreateStaticCondMatrix(
1968  const MatrixKey &mkey)
1969  {
1970  DNekScalBlkMatSharedPtr returnval;
1971 
1972  ASSERTL2(m_metricinfo->GetGtype()
1973  != SpatialDomains::eNoGeomType,
1974  "Geometric information is not set up");
1975 
1976  // set up block matrix system
1977  unsigned int nbdry = NumBndryCoeffs();
1978  unsigned int nint = (unsigned int)(m_ncoeffs - nbdry);
1979  unsigned int exp_size[] = {nbdry,nint};
1980  unsigned int nblks = 2;
1981  returnval = MemoryManager<DNekScalBlkMat>::
1982  AllocateSharedPtr(nblks,nblks,exp_size,exp_size);
1983  //Really need a constructor which takes Arrays
1984  NekDouble factor = 1.0;
1985 
1986  switch (mkey.GetMatrixType())
1987  {
1988  // this can only use stdregions statically condensed system
1989  // for mass matrix
1990  case StdRegions::eMass:
1991  if ((m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1992  ||(mkey.GetNVarCoeff()))
1993  {
1994  factor = 1.0;
1995  goto UseLocRegionsMatrix;
1996  }
1997  else
1998  {
1999  factor = (m_metricinfo->GetJac(GetPointsKeys()))[0];
2000  goto UseStdRegionsMatrix;
2001  }
2002  break;
2003  default: // use Deformed case for both
2004  // regular and deformed geometries
2005  factor = 1.0;
2006  goto UseLocRegionsMatrix;
2007  break;
2008  UseStdRegionsMatrix:
2009  {
2010  NekDouble invfactor = 1.0/factor;
2011  NekDouble one = 1.0;
2012  DNekBlkMatSharedPtr mat = GetStdStaticCondMatrix(mkey);
2013  DNekScalMatSharedPtr Atmp;
2014  DNekMatSharedPtr Asubmat;
2015 
2016  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::
2017  AllocateSharedPtr(factor,Asubmat = mat->GetBlock(0,0)));
2018  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::
2019  AllocateSharedPtr(one,Asubmat = mat->GetBlock(0,1)));
2020  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::
2021  AllocateSharedPtr(factor,Asubmat = mat->GetBlock(1,0)));
2022  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::
2023  AllocateSharedPtr(invfactor,Asubmat = mat->GetBlock(1,1)));
2024  }
2025  break;
2026  UseLocRegionsMatrix:
2027  {
2028  int i,j;
2029  NekDouble invfactor = 1.0/factor;
2030  NekDouble one = 1.0;
2031  DNekScalMat &mat = *GetLocMatrix(mkey);
2032  DNekMatSharedPtr A = MemoryManager<DNekMat>::
2033  AllocateSharedPtr(nbdry,nbdry);
2034  DNekMatSharedPtr B = MemoryManager<DNekMat>::
2035  AllocateSharedPtr(nbdry,nint);
2036  DNekMatSharedPtr C = MemoryManager<DNekMat>::
2037  AllocateSharedPtr(nint,nbdry);
2038  DNekMatSharedPtr D = MemoryManager<DNekMat>::
2039  AllocateSharedPtr(nint,nint);
2040 
2041  Array<OneD,unsigned int> bmap(nbdry);
2042  Array<OneD,unsigned int> imap(nint);
2043  GetBoundaryMap(bmap);
2044  GetInteriorMap(imap);
2045 
2046  for (i = 0; i < nbdry; ++i)
2047  {
2048  for(j = 0; j < nbdry; ++j)
2049  {
2050  (*A)(i,j) = mat(bmap[i],bmap[j]);
2051  }
2052 
2053  for(j = 0; j < nint; ++j)
2054  {
2055  (*B)(i,j) = mat(bmap[i],imap[j]);
2056  }
2057  }
2058 
2059  for (i = 0; i < nint; ++i)
2060  {
2061  for(j = 0; j < nbdry; ++j)
2062  {
2063  (*C)(i,j) = mat(imap[i],bmap[j]);
2064  }
2065 
2066  for(j = 0; j < nint; ++j)
2067  {
2068  (*D)(i,j) = mat(imap[i],imap[j]);
2069  }
2070  }
2071 
2072  // Calculate static condensed system
2073  if(nint)
2074  {
2075  D->Invert();
2076  (*B) = (*B)*(*D);
2077  (*A) = (*A) - (*B)*(*C);
2078  }
2079 
2080  DNekScalMatSharedPtr Atmp;
2081 
2082  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::
2083  AllocateSharedPtr(factor, A));
2084  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::
2085  AllocateSharedPtr(one, B));
2086  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::
2087  AllocateSharedPtr(factor, C));
2088  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::
2089  AllocateSharedPtr(invfactor, D));
2090 
2091  }
2092  }
2093  return returnval;
2094  }
2095 
2096 
2097  DNekScalMatSharedPtr QuadExp::v_GetLocMatrix(const MatrixKey &mkey)
2098  {
2099  return m_matrixManager[mkey];
2100  }
2101 
2102 
2103  DNekScalBlkMatSharedPtr QuadExp::v_GetLocStaticCondMatrix(
2104  const MatrixKey &mkey)
2105  {
2106  return m_staticCondMatrixManager[mkey];
2107  }
2108 
2109  void QuadExp::v_DropLocStaticCondMatrix(const MatrixKey &mkey)
2110  {
2111  m_staticCondMatrixManager.DeleteObject(mkey);
2112  }
2113 
2114 
2115  void QuadExp::v_MassMatrixOp(
2116  const Array<OneD, const NekDouble> &inarray,
2117  Array<OneD,NekDouble> &outarray,
2118  const StdRegions::StdMatrixKey &mkey)
2119  {
2120  StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
2121  }
2122 
2123 
2124  void QuadExp::v_LaplacianMatrixOp(
2125  const Array<OneD, const NekDouble> &inarray,
2126  Array<OneD,NekDouble> &outarray,
2127  const StdRegions::StdMatrixKey &mkey)
2128  {
2129  QuadExp::LaplacianMatrixOp_MatFree(inarray, outarray, mkey);
2130  }
2131 
2132 
2133  void QuadExp::v_LaplacianMatrixOp(
2134  const int k1,
2135  const int k2,
2136  const Array<OneD, const NekDouble> &inarray,
2137  Array<OneD,NekDouble> &outarray,
2138  const StdRegions::StdMatrixKey &mkey)
2139  {
2140  StdExpansion::LaplacianMatrixOp_MatFree(
2141  k1, k2, inarray, outarray, mkey);
2142  }
2143 
2144 
2145  void QuadExp::v_WeakDerivMatrixOp(
2146  const int i,
2147  const Array<OneD, const NekDouble> &inarray,
2148  Array<OneD,NekDouble> &outarray,
2149  const StdRegions::StdMatrixKey &mkey)
2150  {
2151  StdExpansion::WeakDerivMatrixOp_MatFree(i, inarray, outarray, mkey);
2152  }
2153 
2154 
2155  void QuadExp::v_WeakDirectionalDerivMatrixOp(
2156  const Array<OneD, const NekDouble> &inarray,
2157  Array<OneD,NekDouble> &outarray,
2158  const StdRegions::StdMatrixKey &mkey)
2159  {
2160  StdExpansion::WeakDirectionalDerivMatrixOp_MatFree(
2161  inarray, outarray, mkey);
2162  }
2163 
2164 
2165  void QuadExp::v_MassLevelCurvatureMatrixOp(
2166  const Array<OneD, const NekDouble> &inarray,
2167  Array<OneD,NekDouble> &outarray,
2168  const StdRegions::StdMatrixKey &mkey)
2169  {
2170  StdExpansion::MassLevelCurvatureMatrixOp_MatFree(
2171  inarray, outarray, mkey);
2172  }
2173 
2174 
2175  void QuadExp::v_HelmholtzMatrixOp(
2176  const Array<OneD, const NekDouble> &inarray,
2177  Array<OneD,NekDouble> &outarray,
2178  const StdRegions::StdMatrixKey &mkey)
2179  {
2180  QuadExp::HelmholtzMatrixOp_MatFree(inarray, outarray, mkey);
2181  }
2182 
2183 
2184  void QuadExp::v_GeneralMatrixOp_MatOp(
2185  const Array<OneD, const NekDouble> &inarray,
2186  Array<OneD,NekDouble> &outarray,
2187  const StdRegions::StdMatrixKey &mkey)
2188  {
2189  MatrixKey newkey(mkey);
2190  DNekScalMatSharedPtr mat = GetLocMatrix(newkey);
2191 
2192  if (inarray.get() == outarray.get())
2193  {
2194  Array<OneD,NekDouble> tmp(m_ncoeffs);
2195  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
2196 
2197  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs, mat->Scale(),
2198  (mat->GetOwnedMatrix())->GetPtr().get(), m_ncoeffs,
2199  tmp.get(), 1, 0.0, outarray.get(), 1);
2200  }
2201  else
2202  {
2203  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),
2204  (mat->GetOwnedMatrix())->GetPtr().get(), m_ncoeffs,
2205  inarray.get(), 1, 0.0, outarray.get(), 1);
2206  }
2207  }
2208 
2209  void QuadExp::v_ReduceOrderCoeffs(
2210  int numMin,
2211  const Array<OneD, const NekDouble> &inarray,
2212  Array<OneD, NekDouble> &outarray)
2213  {
2214  int n_coeffs = inarray.num_elements();
2215 
2216  Array<OneD, NekDouble> coeff (n_coeffs);
2217  Array<OneD, NekDouble> coeff_tmp(n_coeffs, 0.0);
2218  Array<OneD, NekDouble> tmp, tmp2;
2219 
2220  int nmodes0 = m_base[0]->GetNumModes();
2221  int nmodes1 = m_base[1]->GetNumModes();
2222  int numMax = nmodes0;
2223 
2224  Vmath::Vcopy(n_coeffs,inarray,1,coeff_tmp,1);
2225 
2226  const LibUtilities::PointsKey Pkey0(
2227  nmodes0, LibUtilities::eGaussLobattoLegendre);
2228  const LibUtilities::PointsKey Pkey1(
2229  nmodes1, LibUtilities::eGaussLobattoLegendre);
2230  LibUtilities::BasisKey b0(
2231  m_base[0]->GetBasisType(), nmodes0, Pkey0);
2232  LibUtilities::BasisKey b1(
2233  m_base[1]->GetBasisType(), nmodes1, Pkey1);
2234  LibUtilities::BasisKey bortho0(
2235  LibUtilities::eOrtho_A, nmodes0, Pkey0);
2236  LibUtilities::BasisKey bortho1(
2237  LibUtilities::eOrtho_A, nmodes1, Pkey1);
2238 
2239  LibUtilities::InterpCoeff2D(
2240  b0, b1, coeff_tmp, bortho0, bortho1, coeff);
2241 
2242  Vmath::Zero(n_coeffs, coeff_tmp, 1);
2243 
2244  int cnt = 0;
2245  for (int i = 0; i < numMin+1; ++i)
2246  {
2247  Vmath::Vcopy(numMin,
2248  tmp = coeff+cnt,1,
2249  tmp2 = coeff_tmp+cnt,1);
2250 
2251  cnt = i*numMax;
2252  }
2253 
2254  LibUtilities::InterpCoeff2D(
2255  bortho0, bortho1, coeff_tmp,
2256  b0, b1, outarray);
2257  }
2258 
2259  void QuadExp::v_LaplacianMatrixOp_MatFree_Kernel(
2260  const Array<OneD, const NekDouble> &inarray,
2261  Array<OneD, NekDouble> &outarray,
2262  Array<OneD, NekDouble> &wsp)
2263  {
2264  if (m_metrics.count(eMetricLaplacian00) == 0)
2265  {
2266  ComputeLaplacianMetric();
2267  }
2268 
2269  int nquad0 = m_base[0]->GetNumPoints();
2270  int nquad1 = m_base[1]->GetNumPoints();
2271  int nqtot = nquad0*nquad1;
2272  int nmodes0 = m_base[0]->GetNumModes();
2273  int nmodes1 = m_base[1]->GetNumModes();
2274  int wspsize = max(max(max(nqtot,m_ncoeffs),nquad1*nmodes0),nquad0*nmodes1);
2275 
2276  ASSERTL1(wsp.num_elements() >= 3*wspsize,
2277  "Workspace is of insufficient size.");
2278 
2279  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
2280  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
2281  const Array<OneD, const NekDouble>& dbase0 = m_base[0]->GetDbdata();
2282  const Array<OneD, const NekDouble>& dbase1 = m_base[1]->GetDbdata();
2283  const Array<OneD, const NekDouble>& metric00 = m_metrics[eMetricLaplacian00];
2284  const Array<OneD, const NekDouble>& metric01 = m_metrics[eMetricLaplacian01];
2285  const Array<OneD, const NekDouble>& metric11 = m_metrics[eMetricLaplacian11];
2286 
2287  // Allocate temporary storage
2288  Array<OneD,NekDouble> wsp0(wsp);
2289  Array<OneD,NekDouble> wsp1(wsp+wspsize);
2290  Array<OneD,NekDouble> wsp2(wsp+2*wspsize);
2291 
2292  StdExpansion2D::PhysTensorDeriv(inarray,wsp1,wsp2);
2293 
2294  // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
2295  // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
2296  // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
2297  // especially for this purpose
2298  Vmath::Vvtvvtp(nqtot,&metric00[0],1,&wsp1[0],1,&metric01[0],1,&wsp2[0],1,&wsp0[0],1);
2299  Vmath::Vvtvvtp(nqtot,&metric01[0],1,&wsp1[0],1,&metric11[0],1,&wsp2[0],1,&wsp2[0],1);
2300 
2301  // outarray = m = (D_xi1 * B)^T * k
2302  // wsp1 = n = (D_xi2 * B)^T * l
2303  IProductWRTBase_SumFacKernel(dbase0,base1,wsp0,outarray,wsp1,false,true);
2304  IProductWRTBase_SumFacKernel(base0,dbase1,wsp2,wsp1, wsp0,true,false);
2305 
2306  // outarray = outarray + wsp1
2307  // = L * u_hat
2308  Vmath::Vadd(m_ncoeffs,wsp1.get(),1,outarray.get(),1,outarray.get(),1);
2309  }
2310 
2311  void QuadExp::v_ComputeLaplacianMetric()
2312  {
2313  if (m_metrics.count(eMetricQuadrature) == 0)
2314  {
2315  ComputeQuadratureMetric();
2316  }
2317 
2318  const SpatialDomains::GeomType type = m_metricinfo->GetGtype();
2319  const unsigned int nqtot = GetTotPoints();
2320  const unsigned int dim = 2;
2321  const MetricType m[3][3] = { {eMetricLaplacian00, eMetricLaplacian01, eMetricLaplacian02},
2322  {eMetricLaplacian01, eMetricLaplacian11, eMetricLaplacian12},
2323  {eMetricLaplacian02, eMetricLaplacian12, eMetricLaplacian22}
2324  };
2325 
2326  const Array<TwoD, const NekDouble> gmat =
2327  m_metricinfo->GetGmat(GetPointsKeys());
2328  for (unsigned int i = 0; i < dim; ++i)
2329  {
2330  for (unsigned int j = i; j < dim; ++j)
2331  {
2332  m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
2333  if (type == SpatialDomains::eDeformed)
2334  {
2335  Vmath::Vcopy(nqtot, &gmat[i*dim+j][0], 1,
2336  &m_metrics[m[i][j]][0], 1);
2337  }
2338  else
2339  {
2340  Vmath::Fill(nqtot, gmat[i*dim+j][0],
2341  &m_metrics[m[i][j]][0], 1);
2342  }
2343  MultiplyByQuadratureMetric(m_metrics[m[i][j]],
2344  m_metrics[m[i][j]]);
2345 
2346  }
2347  }
2348  }
2349 
2350  void QuadExp::v_SVVLaplacianFilter(
2351  Array<OneD, NekDouble> &array,
2352  const StdRegions::StdMatrixKey &mkey)
2353  {
2354  int nq = GetTotPoints();
2355 
2356  // Calculate sqrt of the Jacobian
2357  Array<OneD, const NekDouble> jac =
2358  m_metricinfo->GetJac(GetPointsKeys());
2359  Array<OneD, NekDouble> sqrt_jac(nq);
2360  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
2361  {
2362  Vmath::Vsqrt(nq,jac,1,sqrt_jac,1);
2363  }
2364  else
2365  {
2366  Vmath::Fill(nq,sqrt(jac[0]),sqrt_jac,1);
2367  }
2368 
2369  // Multiply array by sqrt(Jac)
2370  Vmath::Vmul(nq,sqrt_jac,1,array,1,array,1);
2371 
2372  // Apply std region filter
2373  StdQuadExp::v_SVVLaplacianFilter( array, mkey);
2374 
2375  // Divide by sqrt(Jac)
2376  Vmath::Vdiv(nq,array,1,sqrt_jac,1,array,1);
2377  }
2378 
2379  }//end of namespace
2380 }//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:2209
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:1603
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:2115
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:1609
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:2109
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:1616
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:1967
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:2097
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:2184
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:2175
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:1597
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:2103
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:1653
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:2124
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:2165
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:2145
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:2259
Vmath.hpp
Nektar::LocalRegions::QuadExp::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: QuadExp.cpp:1621
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:2155
Nektar::LocalRegions::QuadExp::v_ComputeLaplacianMetric
virtual void v_ComputeLaplacianMetric()
Definition: QuadExp.cpp:2311
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:1642
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:2350
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 Wed Mar 8 2017 12:29:00 for Nektar++ by   doxygen 1.8.8