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