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.size())
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.size())
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.size())
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.size())
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.size())
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.size())
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.size())
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] = GetTraceOrient(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  GetTraceToElementMap(i,mapArray,signArray,orient[i]);
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  Array<OneD, NekDouble> tmp1(2*nqtot+m_ncoeffs+nmodes0*nquad1);
502  Array<OneD, NekDouble> tmp2(tmp1 + nqtot);
503  Array<OneD, NekDouble> tmp3(tmp1 + 2*nqtot);
504  Array<OneD, NekDouble> tmp4(tmp1 + 2*nqtot+m_ncoeffs);
505 
506  Array<OneD, Array<OneD, NekDouble>> tmp2D{2};
507  tmp2D[0] = tmp1;
508  tmp2D[1] = tmp2;
509 
510  QuadExp::v_AlignVectorToCollapsedDir(dir, inarray, tmp2D);
511 
512  MultiplyByQuadratureMetric(tmp1,tmp1);
513  MultiplyByQuadratureMetric(tmp2,tmp2);
514 
515  IProductWRTBase_SumFacKernel(
516  m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
517  tmp1, tmp3, tmp4, false, true);
518  IProductWRTBase_SumFacKernel(
519  m_base[0]->GetBdata() , m_base[1]->GetDbdata(),
520  tmp2, outarray, tmp4, true, false);
521  Vmath::Vadd(m_ncoeffs, tmp3, 1, outarray, 1, outarray, 1);
522  }
523 
524  void QuadExp::v_AlignVectorToCollapsedDir(
525  const int dir,
526  const Array<OneD, const NekDouble> &inarray,
527  Array<OneD, Array<OneD, NekDouble> > &outarray)
528  {
529  ASSERTL1((dir==0) || (dir==1) || (dir==2),
530  "Invalid direction.");
531  ASSERTL1((dir==2) ? (m_geom->GetCoordim() ==3):true,
532  "Invalid direction.");
533 
534  int nquad0 = m_base[0]->GetNumPoints();
535  int nquad1 = m_base[1]->GetNumPoints();
536  int nqtot = nquad0*nquad1;
537  int nmodes0 = m_base[0]->GetNumModes();
538 
539  const Array<TwoD, const NekDouble>& df =
540  m_metricinfo->GetDerivFactors(GetPointsKeys());
541 
542  Array<OneD, NekDouble> tmp1 = outarray[0];
543  Array<OneD, NekDouble> tmp2 = outarray[1];
544  Array<OneD, NekDouble> tmp3(m_ncoeffs);
545  Array<OneD, NekDouble> tmp4(nmodes0*nquad1);
546 
547  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
548  {
549  Vmath::Vmul(nqtot,
550  &df[2*dir][0], 1,
551  inarray.get(), 1,
552  tmp1.get(), 1);
553  Vmath::Vmul(nqtot,
554  &df[2*dir+1][0], 1,
555  inarray.get(), 1,
556  tmp2.get(),1);
557  }
558  else
559  {
560  Vmath::Smul(nqtot,
561  df[2*dir][0], inarray.get(), 1,
562  tmp1.get(), 1);
563  Vmath::Smul(nqtot,
564  df[2*dir+1][0], inarray.get(), 1,
565  tmp2.get(), 1);
566  }
567  }
568 
569  void QuadExp::v_IProductWRTDerivBase_MatOp(
570  const int dir,
571  const Array<OneD, const NekDouble>& inarray,
572  Array<OneD, NekDouble> &outarray)
573  {
574  int nq = GetTotPoints();
575  StdRegions::MatrixType mtype = StdRegions::eIProductWRTDerivBase0;
576 
577  switch (dir)
578  {
579  case 0:
580  {
581  mtype = StdRegions::eIProductWRTDerivBase0;
582  }
583  break;
584  case 1:
585  {
586  mtype = StdRegions::eIProductWRTDerivBase1;
587  }
588  break;
589  case 2:
590  {
591  mtype = StdRegions::eIProductWRTDerivBase2;
592  }
593  break;
594  default:
595  {
596  ASSERTL1(false,"input dir is out of range");
597  }
598  break;
599  }
600 
601  MatrixKey iprodmatkey(mtype,DetShapeType(),*this);
602  DNekScalMatSharedPtr iprodmat = m_matrixManager[iprodmatkey];
603 
604  Blas::Dgemv('N', m_ncoeffs, nq, iprodmat->Scale(),
605  (iprodmat->GetOwnedMatrix())->GetPtr().get(),
606  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
607  }
608 
609 
610  void QuadExp::v_NormVectorIProductWRTBase(
611  const Array<OneD, const NekDouble> &Fx,
612  const Array<OneD, const NekDouble> &Fy,
613  const Array<OneD, const NekDouble> &Fz,
614  Array<OneD, NekDouble> &outarray)
615  {
616  int nq = m_base[0]->GetNumPoints()*m_base[1]->GetNumPoints();
617  Array<OneD, NekDouble> Fn(nq);
618 
619  const Array<OneD, const Array<OneD, NekDouble> > &normals =
620  GetLeftAdjacentElementExp()->GetTraceNormal(
621  GetLeftAdjacentElementTrace());
622 
623  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
624  {
625  Vmath::Vvtvvtp(nq,&normals[0][0],1,&Fx[0],1,
626  &normals[1][0],1,&Fy[0],1,&Fn[0],1);
627  Vmath::Vvtvp (nq,&normals[2][0],1,&Fz[0],1,&Fn[0],1,&Fn[0],1);
628  }
629  else
630  {
631  Vmath::Svtsvtp(nq,normals[0][0],&Fx[0],1,
632  normals[1][0],&Fy[0],1,&Fn[0],1);
633  Vmath::Svtvp (nq,normals[2][0],&Fz[0],1,&Fn[0],1,&Fn[0],1);
634  }
635 
636  IProductWRTBase(Fn,outarray);
637  }
638 
639  void QuadExp::v_NormVectorIProductWRTBase(
640  const Array<OneD, const Array<OneD, NekDouble> > &Fvec,
641  Array<OneD, NekDouble> &outarray)
642  {
643  NormVectorIProductWRTBase(Fvec[0], Fvec[1], Fvec[2], outarray);
644  }
645 
646  StdRegions::StdExpansionSharedPtr QuadExp::v_GetStdExp(void) const
647  {
648  return MemoryManager<StdRegions::StdQuadExp>
649  ::AllocateSharedPtr(m_base[0]->GetBasisKey(),
650  m_base[1]->GetBasisKey());
651  }
652 
653 
654  StdRegions::StdExpansionSharedPtr QuadExp::v_GetLinStdExp(void) const
655  {
656  LibUtilities::BasisKey bkey0(m_base[0]->GetBasisType(),
657  2, m_base[0]->GetPointsKey());
658  LibUtilities::BasisKey bkey1(m_base[1]->GetBasisType(),
659  2, m_base[1]->GetPointsKey());
660 
661  return MemoryManager<StdRegions::StdQuadExp>
662  ::AllocateSharedPtr( bkey0, bkey1);
663  }
664 
665  void QuadExp::v_GetCoords(
666  Array<OneD, NekDouble> &coords_0,
667  Array<OneD, NekDouble> &coords_1,
668  Array<OneD, NekDouble> &coords_2)
669  {
670  Expansion::v_GetCoords(coords_0, coords_1, coords_2);
671  }
672 
673 
674  void QuadExp::v_GetCoord(const Array<OneD, const NekDouble> &Lcoords,
675  Array<OneD,NekDouble> &coords)
676  {
677  int i;
678 
679  ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[1] <= 1.0 &&
680  Lcoords[1] >= -1.0 && Lcoords[1] <=1.0,
681  "Local coordinates are not in region [-1,1]");
682 
683  m_geom->FillGeom();
684  for (i = 0; i < m_geom->GetCoordim(); ++i)
685  {
686  coords[i] = m_geom->GetCoord(i,Lcoords);
687  }
688  }
689 
690 
691 
692  /**
693  * Given the local cartesian coordinate \a Lcoord evaluate the
694  * value of physvals at this point by calling through to the
695  * StdExpansion method
696  */
697  NekDouble QuadExp::v_StdPhysEvaluate(
698  const Array<OneD, const NekDouble> &Lcoord,
699  const Array<OneD, const NekDouble> &physvals)
700  {
701  // Evaluate point in local (eta) coordinates.
702  return StdQuadExp::v_PhysEvaluate(Lcoord,physvals);
703  }
704 
705  NekDouble QuadExp::v_PhysEvaluate(
706  const Array<OneD, const NekDouble> &coord,
707  const Array<OneD, const NekDouble> &physvals)
708  {
709  Array<OneD,NekDouble> Lcoord = Array<OneD, NekDouble>(2);
710 
711  ASSERTL0(m_geom,"m_geom not defined");
712  m_geom->GetLocCoords(coord,Lcoord);
713 
714  return StdQuadExp::v_PhysEvaluate(Lcoord, physvals);
715  }
716 
717 
718  // Get edge values from the 2D Phys space along an edge
719  // following a counter clockwise edge convention for definition
720  // of edgedir, Note that point distribution is given by QuadExp.
721  void QuadExp::v_GetEdgePhysVals(
722  const int edge,
723  const Array<OneD, const NekDouble> &inarray,
724  Array<OneD,NekDouble> &outarray)
725  {
726  int nquad0 = m_base[0]->GetNumPoints();
727  int nquad1 = m_base[1]->GetNumPoints();
728 
729  StdRegions::Orientation edgedir = GetTraceOrient(edge);
730  switch(edge)
731  {
732  case 0:
733  if (edgedir == StdRegions::eForwards)
734  {
735  Vmath::Vcopy(nquad0,&(inarray[0]),1,&(outarray[0]),1);
736  }
737  else
738  {
739  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0-1),-1,
740  &(outarray[0]),1);
741  }
742  break;
743  case 1:
744  if (edgedir == StdRegions::eForwards)
745  {
746  Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0-1),nquad0,
747  &(outarray[0]),1);
748  }
749  else
750  {
751  Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0*nquad1-1),
752  -nquad0, &(outarray[0]),1);
753  }
754  break;
755  case 2:
756  if (edgedir == StdRegions::eForwards)
757  {
758  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1-1),-1,
759  &(outarray[0]),1);
760  }
761  else
762  {
763  Vmath::Vcopy(nquad0,&(inarray[0])+nquad0*(nquad1-1),1,
764  &(outarray[0]),1);
765  }
766  break;
767  case 3:
768  if (edgedir == StdRegions::eForwards)
769  {
770  Vmath::Vcopy(nquad1,&(inarray[0]) + nquad0*(nquad1-1),
771  -nquad0,&(outarray[0]),1);
772  }
773  else
774  {
775  Vmath::Vcopy(nquad1,&(inarray[0]),nquad0,
776  &(outarray[0]),1);
777  }
778  break;
779  default:
780  ASSERTL0(false,"edge value (< 3) is out of range");
781  break;
782  }
783  }
784 
785 
786  void QuadExp::v_GetTracePhysVals(
787  const int edge,
788  const StdRegions::StdExpansionSharedPtr &EdgeExp,
789  const Array<OneD, const NekDouble> &inarray,
790  Array<OneD,NekDouble> &outarray,
791  StdRegions::Orientation orient)
792  {
793  int nquad0 = m_base[0]->GetNumPoints();
794  int nquad1 = m_base[1]->GetNumPoints();
795 
796  // Implementation for all the basis except Gauss points
797  if (m_base[0]->GetPointsType() !=
798  LibUtilities::eGaussGaussLegendre &&
799  m_base[1]->GetPointsType() !=
800  LibUtilities::eGaussGaussLegendre)
801  {
802  switch (edge)
803  {
804  case 0:
805  Vmath::Vcopy(nquad0,&(inarray[0]),1,&(outarray[0]),1);
806  break;
807  case 1:
808  Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0-1),
809  nquad0,&(outarray[0]),1);
810  break;
811  case 2:
812  Vmath::Vcopy(nquad0,&(inarray[0])+nquad0*(nquad1-1),1,
813  &(outarray[0]),1);
814  break;
815  case 3:
816  Vmath::Vcopy(nquad1,&(inarray[0]),nquad0,
817  &(outarray[0]),1);
818  break;
819  default:
820  ASSERTL0(false,"edge value (< 3) is out of range");
821  break;
822  }
823  }
824  else
825  {
826  QuadExp::v_GetEdgeInterpVals(edge, inarray, outarray);
827  }
828 
829  // Interpolate if required
830  if (m_base[edge%2]->GetPointsKey() !=
831  EdgeExp->GetBasis(0)->GetPointsKey())
832  {
833  Array<OneD,NekDouble> outtmp(max(nquad0,nquad1));
834 
835  outtmp = outarray;
836 
837  LibUtilities::Interp1D(
838  m_base[edge%2]->GetPointsKey(), outtmp,
839  EdgeExp->GetBasis(0)->GetPointsKey(), outarray);
840  }
841 
842  if (orient == StdRegions::eNoOrientation)
843  {
844  orient = GetTraceOrient(edge);
845  }
846  //Reverse data if necessary
847  if(orient == StdRegions::eBackwards)
848  {
849  Vmath::Reverse(EdgeExp->GetNumPoints(0),&outarray[0], 1,
850  &outarray[0], 1);
851  }
852  }
853 
854  void QuadExp::v_GetEdgeInterpVals(const int edge,
855  const Array<OneD, const NekDouble> &inarray,
856  Array<OneD,NekDouble> &outarray)
857  {
858  int i;
859  int nq0 = m_base[0]->GetNumPoints();
860  int nq1 = m_base[1]->GetNumPoints();
861 
862  StdRegions::ConstFactorMap factors;
863  factors[StdRegions::eFactorGaussEdge] = edge;
864 
865  StdRegions::StdMatrixKey key(
866  StdRegions::eInterpGauss,
867  DetShapeType(),*this,factors);
868 
869  DNekScalMatSharedPtr mat_gauss = m_matrixManager[key];
870 
871  switch (edge)
872  {
873  case 0:
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 1:
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  case 2:
894  {
895  for (i = 0; i < nq0; i++)
896  {
897  outarray[i] = Blas::Ddot(
898  nq1, mat_gauss->GetOwnedMatrix()->GetPtr().get(),
899  1, &inarray[i], nq0);
900  }
901  break;
902  }
903  case 3:
904  {
905  for (i = 0; i < nq1; i++)
906  {
907  outarray[i] = Blas::Ddot(
908  nq0, mat_gauss->GetOwnedMatrix()->GetPtr().get(),
909  1, &inarray[i * nq0], 1);
910  }
911  break;
912  }
913  default:
914  ASSERTL0(false, "edge value (< 3) is out of range");
915  break;
916  }
917  }
918 
919 
920  void QuadExp::v_GetTracePhysMap(
921  const int edge,
922  Array<OneD, int> &outarray)
923  {
924  int nquad0 = m_base[0]->GetNumPoints();
925  int nquad1 = m_base[1]->GetNumPoints();
926 
927  // Get points in Cartesian orientation
928  switch (edge)
929  {
930  case 0:
931  outarray = Array<OneD, int>(nquad0);
932  for (int i = 0; i < nquad0; ++i)
933  {
934  outarray[i] = i;
935  }
936  break;
937  case 1:
938  outarray = Array<OneD, int>(nquad1);
939  for (int i = 0; i < nquad1; ++i)
940  {
941  outarray[i] = (nquad0-1) + i*nquad0;
942  }
943  break;
944  case 2:
945  outarray = Array<OneD, int>(nquad0);
946  for (int i = 0; i < nquad0; ++i)
947  {
948  outarray[i] = i + nquad0*(nquad1-1);
949  }
950  break;
951  case 3:
952  outarray = Array<OneD, int>(nquad1);
953  for (int i = 0; i < nquad1; ++i)
954  {
955  outarray[i] = i*nquad0;
956  }
957  break;
958  default:
959  ASSERTL0(false, "edge value (< 3) is out of range");
960  break;
961  }
962 
963  }
964 
965 
966 
967 
968  void QuadExp::v_GetTraceQFactors(
969  const int edge,
970  Array<OneD, NekDouble> &outarray)
971  {
972  int i;
973  int nquad0 = m_base[0]->GetNumPoints();
974  int nquad1 = m_base[1]->GetNumPoints();
975 
976  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
977  const Array<OneD, const NekDouble>& jac = m_metricinfo->GetJac(ptsKeys);
978  const Array<TwoD, const NekDouble>& df = m_metricinfo->GetDerivFactors(ptsKeys);
979 
980  Array<OneD, NekDouble> j (max(nquad0, nquad1), 0.0);
981  Array<OneD, NekDouble> g0(max(nquad0, nquad1), 0.0);
982  Array<OneD, NekDouble> g1(max(nquad0, nquad1), 0.0);
983  Array<OneD, NekDouble> g2(max(nquad0, nquad1), 0.0);
984  Array<OneD, NekDouble> g3(max(nquad0, nquad1), 0.0);
985 
986  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
987  {
988  // Implementation for all the basis except Gauss points
989  if (m_base[0]->GetPointsType()
990  != LibUtilities::eGaussGaussLegendre
991  && m_base[1]->GetPointsType() !=
992  LibUtilities::eGaussGaussLegendre)
993  {
994  switch (edge)
995  {
996  case 0:
997  Vmath::Vcopy(nquad0, &(df[1][0]),
998  1, &(g1[0]), 1);
999  Vmath::Vcopy(nquad0, &(df[3][0]),
1000  1, &(g3[0]), 1);
1001  Vmath::Vcopy(nquad0, &(jac[0]),1, &(j[0]), 1);
1002 
1003  for (i = 0; i < nquad0; ++i)
1004  {
1005  outarray[i] = j[i]*sqrt(g1[i]*g1[i]
1006  + g3[i]*g3[i]);
1007  }
1008  break;
1009  case 1:
1010  Vmath::Vcopy(nquad1,
1011  &(df[0][0])+(nquad0-1), nquad0,
1012  &(g0[0]), 1);
1013 
1014  Vmath::Vcopy(nquad1,
1015  &(df[2][0])+(nquad0-1), nquad0,
1016  &(g2[0]), 1);
1017 
1018  Vmath::Vcopy(nquad1,
1019  &(jac[0])+(nquad0-1), nquad0,
1020  &(j[0]), 1);
1021 
1022  for (i = 0; i < nquad1; ++i)
1023  {
1024  outarray[i] = j[i]*sqrt(g0[i]*g0[i] +
1025  g2[i]*g2[i]);
1026  }
1027  break;
1028  case 2:
1029 
1030  Vmath::Vcopy(nquad0,
1031  &(df[1][0])+(nquad0*(nquad1-1)), 1,
1032  &(g1[0]), 1);
1033 
1034  Vmath::Vcopy(nquad0,
1035  &(df[3][0])+(nquad0*(nquad1-1)), 1,
1036  &(g3[0]), 1);
1037 
1038  Vmath::Vcopy(nquad0,
1039  &(jac[0])+(nquad0*(nquad1-1)), 1,
1040  &(j[0]), 1);
1041 
1042  for (i = 0; i < nquad0; ++i)
1043  {
1044  outarray[i] =
1045  j[i]*sqrt(g1[i]*g1[i]+ g3[i]*g3[i]);
1046  }
1047  break;
1048  case 3:
1049 
1050  Vmath::Vcopy(nquad1, &(df[0][0]), nquad0,&(g0[0]), 1);
1051  Vmath::Vcopy(nquad1, &(df[2][0]), nquad0,&(g2[0]), 1);
1052  Vmath::Vcopy(nquad1, &(jac[0]), nquad0, &(j[0]), 1);
1053 
1054  for (i = 0; i < nquad1; ++i)
1055  {
1056  outarray[i] = j[i]*sqrt(g0[i]*g0[i] +
1057  g2[i]*g2[i]);
1058  }
1059  break;
1060  default:
1061  ASSERTL0(false,"edge value (< 3) is out of range");
1062  break;
1063  }
1064  }
1065  else
1066  {
1067  int nqtot = nquad0 * nquad1;
1068  Array<OneD, NekDouble> tmp_gmat0(nqtot, 0.0);
1069  Array<OneD, NekDouble> tmp_gmat1(nqtot, 0.0);
1070  Array<OneD, NekDouble> tmp_gmat2(nqtot, 0.0);
1071  Array<OneD, NekDouble> tmp_gmat3(nqtot, 0.0);
1072  Array<OneD, NekDouble> g0_edge(max(nquad0, nquad1), 0.0);
1073  Array<OneD, NekDouble> g1_edge(max(nquad0, nquad1), 0.0);
1074  Array<OneD, NekDouble> g2_edge(max(nquad0, nquad1), 0.0);
1075  Array<OneD, NekDouble> g3_edge(max(nquad0, nquad1), 0.0);
1076  Array<OneD, NekDouble> jac_edge(max(nquad0, nquad1), 0.0);
1077 
1078  switch (edge)
1079  {
1080  case 0:
1081  Vmath::Vmul(nqtot,&(df[1][0]),1,&jac[0],1,
1082  &(tmp_gmat1[0]),1);
1083  Vmath::Vmul(nqtot,&(df[3][0]),1,&jac[0],1,
1084  &(tmp_gmat3[0]),1);
1085  QuadExp::v_GetEdgeInterpVals(
1086  edge, tmp_gmat1, g1_edge);
1087  QuadExp::v_GetEdgeInterpVals(
1088  edge, tmp_gmat3, g3_edge);
1089 
1090  for (i = 0; i < nquad0; ++i)
1091  {
1092  outarray[i] = sqrt(g1_edge[i]*g1_edge[i] +
1093  g3_edge[i]*g3_edge[i]);
1094  }
1095  break;
1096 
1097  case 1:
1098  Vmath::Vmul(nqtot,
1099  &(df[0][0]), 1,
1100  &jac[0], 1,
1101  &(tmp_gmat0[0]), 1);
1102  Vmath::Vmul(nqtot,
1103  &(df[2][0]), 1,
1104  &jac[0], 1,
1105  &(tmp_gmat2[0]),
1106  1);
1107  QuadExp::v_GetEdgeInterpVals(
1108  edge, tmp_gmat0, g0_edge);
1109  QuadExp::v_GetEdgeInterpVals(
1110  edge, tmp_gmat2, g2_edge);
1111 
1112  for (i = 0; i < nquad1; ++i)
1113  {
1114  outarray[i] = sqrt(g0_edge[i]*g0_edge[i]
1115  + g2_edge[i]*g2_edge[i]);
1116  }
1117 
1118  break;
1119  case 2:
1120 
1121  Vmath::Vmul(nqtot,
1122  &(df[1][0]), 1,
1123  &jac[0], 1,
1124  &(tmp_gmat1[0]), 1);
1125  Vmath::Vmul(nqtot,
1126  &(df[3][0]), 1,
1127  &jac[0], 1,
1128  &(tmp_gmat3[0]),1);
1129  QuadExp::v_GetEdgeInterpVals(
1130  edge, tmp_gmat1, g1_edge);
1131  QuadExp::v_GetEdgeInterpVals(
1132  edge, tmp_gmat3, g3_edge);
1133 
1134 
1135  for (i = 0; i < nquad0; ++i)
1136  {
1137  outarray[i] = sqrt(g1_edge[i]*g1_edge[i]
1138  + g3_edge[i]*g3_edge[i]);
1139  }
1140 
1141  Vmath::Reverse(nquad0,&outarray[0],1,&outarray[0],1);
1142 
1143  break;
1144  case 3:
1145  Vmath::Vmul(nqtot,
1146  &(df[0][0]), 1,
1147  &jac[0], 1,
1148  &(tmp_gmat0[0]), 1);
1149  Vmath::Vmul(nqtot,
1150  &(df[2][0]),1,
1151  &jac[0], 1,
1152  &(tmp_gmat2[0]),1);
1153  QuadExp::v_GetEdgeInterpVals(
1154  edge, tmp_gmat0, g0_edge);
1155  QuadExp::v_GetEdgeInterpVals(
1156  edge, tmp_gmat2, g2_edge);
1157 
1158 
1159  for (i = 0; i < nquad1; ++i)
1160  {
1161  outarray[i] = sqrt(g0_edge[i]*g0_edge[i] +
1162  g2_edge[i]*g2_edge[i]);
1163  }
1164 
1165  Vmath::Reverse(nquad1,
1166  &outarray[0], 1,
1167  &outarray[0], 1);
1168 
1169  break;
1170  default:
1171  ASSERTL0(false,"edge value (< 3) is out of range");
1172  break;
1173  }
1174  }
1175  }
1176  else
1177  {
1178 
1179  switch (edge)
1180  {
1181  case 0:
1182 
1183 
1184 
1185  for (i = 0; i < nquad0; ++i)
1186  {
1187  outarray[i] = jac[0]*sqrt(df[1][0]*df[1][0] +
1188  df[3][0]*df[3][0]);
1189  }
1190  break;
1191  case 1:
1192  for (i = 0; i < nquad1; ++i)
1193  {
1194  outarray[i] = jac[0]*sqrt(df[0][0]*df[0][0] +
1195  df[2][0]*df[2][0]);
1196  }
1197  break;
1198  case 2:
1199  for (i = 0; i < nquad0; ++i)
1200  {
1201  outarray[i] = jac[0]*sqrt(df[1][0]*df[1][0] +
1202  df[3][0]*df[3][0]);
1203  }
1204  break;
1205  case 3:
1206  for (i = 0; i < nquad1; ++i)
1207  {
1208  outarray[i] = jac[0]*sqrt(df[0][0]*df[0][0] +
1209  df[2][0]*df[2][0]);
1210  }
1211  break;
1212  default:
1213  ASSERTL0(false,"edge value (< 3) is out of range");
1214  break;
1215  }
1216  }
1217  }
1218 
1219 
1220  void QuadExp::v_ComputeTraceNormal(const int edge)
1221  {
1222  int i;
1223  const SpatialDomains::GeomFactorsSharedPtr & geomFactors =
1224  GetGeom()->GetMetricInfo();
1225  SpatialDomains::GeomType type = geomFactors->GetGtype();
1226 
1227  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
1228  for(i = 0; i < ptsKeys.size(); ++i)
1229  {
1230  // Need at least 2 points for computing normals
1231  if (ptsKeys[i].GetNumPoints() == 1)
1232  {
1233  LibUtilities::PointsKey pKey(2, ptsKeys[i].GetPointsType());
1234  ptsKeys[i] = pKey;
1235  }
1236  }
1237 
1238  const Array<TwoD, const NekDouble> & df = geomFactors->GetDerivFactors(ptsKeys);
1239  const Array<OneD, const NekDouble> & jac = geomFactors->GetJac(ptsKeys);
1240  int nqe;
1241  if (edge == 0 || edge == 2)
1242  {
1243  nqe = m_base[0]->GetNumPoints();
1244  }
1245  else
1246  {
1247  nqe = m_base[1]->GetNumPoints();
1248  }
1249  int vCoordDim = GetCoordim();
1250 
1251  m_edgeNormals[edge] = Array<OneD, Array<OneD, NekDouble> >
1252  (vCoordDim);
1253  Array<OneD, Array<OneD, NekDouble> > &normal = m_edgeNormals[edge];
1254  for (i = 0; i < vCoordDim; ++i)
1255  {
1256  normal[i] = Array<OneD, NekDouble>(nqe);
1257  }
1258 
1259  size_t nqb = nqe;
1260  size_t nbnd= edge;
1261  m_elmtBndNormDirElmtLen[nbnd] = Array<OneD, NekDouble> {nqb, 0.0};
1262  Array<OneD, NekDouble> &length = m_elmtBndNormDirElmtLen[nbnd];
1263 
1264  // Regular geometry case
1265  if ((type == SpatialDomains::eRegular)||
1266  (type == SpatialDomains::eMovingRegular))
1267  {
1268  NekDouble fac;
1269  // Set up normals
1270  switch (edge)
1271  {
1272  case 0:
1273  for (i = 0; i < vCoordDim; ++i)
1274  {
1275  Vmath::Fill(nqe, -df[2*i+1][0], normal[i], 1);
1276  }
1277  break;
1278  case 1:
1279  for (i = 0; i < vCoordDim; ++i)
1280  {
1281  Vmath::Fill(nqe, df[2*i][0], normal[i], 1);
1282  }
1283  break;
1284  case 2:
1285  for (i = 0; i < vCoordDim; ++i)
1286  {
1287  Vmath::Fill(nqe, df[2*i+1][0], normal[i], 1);
1288  }
1289  break;
1290  case 3:
1291  for (i = 0; i < vCoordDim; ++i)
1292  {
1293  Vmath::Fill(nqe, -df[2*i][0], normal[i], 1);
1294  }
1295  break;
1296  default:
1297  ASSERTL0(false, "edge is out of range (edge < 4)");
1298  }
1299 
1300  // normalise
1301  fac = 0.0;
1302  for (i =0 ; i < vCoordDim; ++i)
1303  {
1304  fac += normal[i][0]*normal[i][0];
1305  }
1306  fac = 1.0/sqrt(fac);
1307 
1308  Vmath::Fill(nqb, fac, length, 1);
1309 
1310  for (i = 0; i < vCoordDim; ++i)
1311  {
1312  Vmath::Smul(nqe, fac, normal[i], 1,normal[i], 1);
1313  }
1314  }
1315  else // Set up deformed normals
1316  {
1317  int j;
1318 
1319  int nquad0 = ptsKeys[0].GetNumPoints();
1320  int nquad1 = ptsKeys[1].GetNumPoints();
1321 
1322  LibUtilities::PointsKey from_key;
1323 
1324  Array<OneD,NekDouble> normals(vCoordDim*max(nquad0,nquad1),0.0);
1325  Array<OneD,NekDouble> edgejac(vCoordDim*max(nquad0,nquad1),0.0);
1326 
1327  // Extract Jacobian along edges and recover local
1328  // derivates (dx/dr) for polynomial interpolation by
1329  // multiplying m_gmat by jacobian
1330 
1331  // Implementation for all the basis except Gauss points
1332  if (m_base[0]->GetPointsType() !=
1333  LibUtilities::eGaussGaussLegendre
1334  && m_base[1]->GetPointsType() !=
1335  LibUtilities::eGaussGaussLegendre)
1336  {
1337  switch (edge)
1338  {
1339  case 0:
1340  for (j = 0; j < nquad0; ++j)
1341  {
1342  edgejac[j] = jac[j];
1343  for (i = 0; i < vCoordDim; ++i)
1344  {
1345  normals[i*nquad0+j] =
1346  -df[2*i+1][j]*edgejac[j];
1347  }
1348  }
1349  from_key = ptsKeys[0];
1350  break;
1351  case 1:
1352  for (j = 0; j < nquad1; ++j)
1353  {
1354  edgejac[j] = jac[nquad0*j+nquad0-1];
1355  for (i = 0; i < vCoordDim; ++i)
1356  {
1357  normals[i*nquad1+j] =
1358  df[2*i][nquad0*j + nquad0-1]
1359  *edgejac[j];
1360  }
1361  }
1362  from_key = ptsKeys[1];
1363  break;
1364  case 2:
1365  for (j = 0; j < nquad0; ++j)
1366  {
1367  edgejac[j] = jac[nquad0*(nquad1-1)+j];
1368  for (i = 0; i < vCoordDim; ++i)
1369  {
1370  normals[i*nquad0+j] =
1371  (df[2*i+1][nquad0*(nquad1-1)+j])
1372  *edgejac[j];
1373  }
1374  }
1375  from_key = ptsKeys[0];
1376  break;
1377  case 3:
1378  for (j = 0; j < nquad1; ++j)
1379  {
1380  edgejac[j] = jac[nquad0*j];
1381  for (i = 0; i < vCoordDim; ++i)
1382  {
1383  normals[i*nquad1+j] =
1384  -df[2*i][nquad0*j]*edgejac[j];
1385  }
1386  }
1387  from_key = ptsKeys[1];
1388  break;
1389  default:
1390  ASSERTL0(false,"edge is out of range (edge < 3)");
1391  }
1392  }
1393  else
1394  {
1395  int nqtot = nquad0 * nquad1;
1396  Array<OneD, NekDouble> tmp_gmat(nqtot, 0.0);
1397  Array<OneD, NekDouble> tmp_gmat_edge(nqe, 0.0);
1398 
1399  switch (edge)
1400  {
1401  case 0:
1402  for (j = 0; j < nquad0; ++j)
1403  {
1404  for (i = 0; i < vCoordDim; ++i)
1405  {
1406  Vmath::Vmul(nqtot,
1407  &(df[2*i+1][0]), 1,
1408  &jac[0], 1,
1409  &(tmp_gmat[0]), 1);
1410  QuadExp::v_GetEdgeInterpVals(
1411  edge, tmp_gmat, tmp_gmat_edge);
1412  normals[i*nquad0+j] = -tmp_gmat_edge[j];
1413  }
1414  }
1415  from_key = ptsKeys[0];
1416  break;
1417  case 1:
1418  for (j = 0; j < nquad1; ++j)
1419  {
1420  for (i = 0; i < vCoordDim; ++i)
1421  {
1422  Vmath::Vmul(nqtot,
1423  &(df[2*i][0]), 1,
1424  &jac[0], 1,
1425  &(tmp_gmat[0]), 1);
1426  QuadExp::v_GetEdgeInterpVals(
1427  edge, tmp_gmat, tmp_gmat_edge);
1428  normals[i*nquad1+j] = tmp_gmat_edge[j];
1429  }
1430  }
1431  from_key = ptsKeys[1];
1432  break;
1433  case 2:
1434  for (j = 0; j < nquad0; ++j)
1435  {
1436  for (i = 0; i < vCoordDim; ++i)
1437  {
1438  Vmath::Vmul(nqtot,
1439  &(df[2*i+1][0]), 1,
1440  &jac[0], 1,
1441  &(tmp_gmat[0]), 1);
1442  QuadExp::v_GetEdgeInterpVals(
1443  edge, tmp_gmat, tmp_gmat_edge);
1444  normals[i*nquad0+j] = tmp_gmat_edge[j];
1445  }
1446  }
1447  from_key = ptsKeys[0];
1448  break;
1449  case 3:
1450  for (j = 0; j < nquad1; ++j)
1451  {
1452  for (i = 0; i < vCoordDim; ++i)
1453  {
1454  Vmath::Vmul(nqtot,
1455  &(df[2*i][0]), 1,
1456  &jac[0], 1,
1457  &(tmp_gmat[0]) ,1);
1458  QuadExp::v_GetEdgeInterpVals(
1459  edge, tmp_gmat, tmp_gmat_edge);
1460  normals[i*nquad1+j] = -tmp_gmat_edge[j];
1461  }
1462  }
1463  from_key = ptsKeys[1];
1464  break;
1465  default:
1466  ASSERTL0(false,"edge is out of range (edge < 3)");
1467  }
1468  }
1469 
1470  int nq = from_key.GetNumPoints();
1471  Array<OneD,NekDouble> work(nqe,0.0);
1472 
1473  // interpolate Jacobian and invert
1474  LibUtilities::Interp1D(
1475  from_key,jac, m_base[0]->GetPointsKey(), work);
1476  Vmath::Sdiv(nqe,1.0,&work[0],1,&work[0],1);
1477 
1478  // interpolate
1479  for (i = 0; i < GetCoordim(); ++i)
1480  {
1481  LibUtilities::Interp1D(
1482  from_key,&normals[i*nq],
1483  m_base[0]->GetPointsKey(),
1484  &normal[i][0]);
1485  Vmath::Vmul(nqe, work, 1, normal[i], 1, normal[i], 1);
1486  }
1487 
1488  //normalise normal vectors
1489  Vmath::Zero(nqe,work,1);
1490  for (i = 0; i < GetCoordim(); ++i)
1491  {
1492  Vmath::Vvtvp(nqe,
1493  normal[i], 1,
1494  normal[i],1 ,
1495  work, 1,
1496  work, 1);
1497  }
1498 
1499  Vmath::Vsqrt(nqe,work,1,work,1);
1500  Vmath::Sdiv(nqe,1.0,work,1,work,1);
1501 
1502  Vmath::Vcopy(nqb, work, 1, length, 1);
1503 
1504  for (i = 0; i < GetCoordim(); ++i)
1505  {
1506  Vmath::Vmul(nqe, normal[i], 1, work, 1, normal[i], 1);
1507  }
1508  }
1509  if (GetGeom()->GetEorient(edge) == StdRegions::eBackwards)
1510  {
1511  for (i = 0; i < vCoordDim; ++i)
1512  {
1513  if (geomFactors->GetGtype() == SpatialDomains::eDeformed)
1514  {
1515  Vmath::Reverse(nqe, normal[i], 1, normal[i],1);
1516  }
1517  }
1518  }
1519  }
1520 
1521  const SpatialDomains::GeomFactorsSharedPtr& QuadExp::v_GetMetricInfo() const
1522  {
1523  return m_metricinfo;
1524  }
1525 
1526 
1527  int QuadExp::v_GetCoordim()
1528  {
1529  return m_geom->GetCoordim();
1530  }
1531 
1532 
1533  void QuadExp::v_ExtractDataToCoeffs(
1534  const NekDouble *data,
1535  const std::vector<unsigned int > &nummodes,
1536  int mode_offset,
1537  NekDouble *coeffs,
1538  std::vector<LibUtilities::BasisType> &fromType)
1539  {
1540  int data_order0 = nummodes[mode_offset];
1541  int fillorder0 = std::min(m_base[0]->GetNumModes(),data_order0);
1542 
1543  int data_order1 = nummodes[mode_offset + 1];
1544  int order1 = m_base[1]->GetNumModes();
1545  int fillorder1 = min(order1,data_order1);
1546 
1547  // Check if same basis
1548  if (fromType[0] != m_base[0]->GetBasisType() ||
1549  fromType[1] != m_base[1]->GetBasisType())
1550  {
1551  // Construct a quad with the appropriate basis type at our
1552  // quadrature points, and one more to do a forwards
1553  // transform. We can then copy the output to coeffs.
1554  StdRegions::StdQuadExp tmpQuad(
1555  LibUtilities::BasisKey(
1556  fromType[0], data_order0, m_base[0]->GetPointsKey()),
1557  LibUtilities::BasisKey(
1558  fromType[1], data_order1, m_base[1]->GetPointsKey()));
1559  StdRegions::StdQuadExp tmpQuad2(m_base[0]->GetBasisKey(),
1560  m_base[1]->GetBasisKey());
1561 
1562  Array<OneD, const NekDouble> tmpData(tmpQuad.GetNcoeffs(), data);
1563  Array<OneD, NekDouble> tmpBwd(tmpQuad2.GetTotPoints());
1564  Array<OneD, NekDouble> tmpOut(tmpQuad2.GetNcoeffs());
1565 
1566  tmpQuad.BwdTrans(tmpData, tmpBwd);
1567  tmpQuad2.FwdTrans(tmpBwd, tmpOut);
1568  Vmath::Vcopy(tmpOut.size(), &tmpOut[0], 1, coeffs, 1);
1569 
1570  return;
1571  }
1572 
1573  switch (m_base[0]->GetBasisType())
1574  {
1575  case LibUtilities::eModified_A:
1576  {
1577  int i;
1578  int cnt = 0;
1579  int cnt1 = 0;
1580 
1581  ASSERTL1(m_base[1]->GetBasisType() ==
1582  LibUtilities::eModified_A,
1583  "Extraction routine not set up for this basis");
1584 
1585  Vmath::Zero(m_ncoeffs,coeffs,1);
1586  for (i = 0; i < fillorder0; ++i)
1587  {
1588  Vmath::Vcopy(fillorder1, data + cnt, 1, coeffs +cnt1, 1);
1589  cnt += data_order1;
1590  cnt1 += order1;
1591  }
1592  }
1593  break;
1594  case LibUtilities::eGLL_Lagrange:
1595  {
1596  LibUtilities::PointsKey
1597  p0(nummodes[0], LibUtilities::eGaussLobattoLegendre);
1598  LibUtilities::PointsKey
1599  p1(nummodes[1], LibUtilities::eGaussLobattoLegendre);
1600  LibUtilities::PointsKey t0(
1601  m_base[0]->GetNumModes(),
1602  LibUtilities::eGaussLobattoLegendre);
1603  LibUtilities::PointsKey t1(
1604  m_base[1]->GetNumModes(),
1605  LibUtilities::eGaussLobattoLegendre);
1606  LibUtilities::Interp2D(p0, p1, data, t0, t1, coeffs);
1607  }
1608  break;
1609  case LibUtilities::eGauss_Lagrange:
1610  {
1611  // Assume that input is also Gll_Lagrange but no way to check;
1612  LibUtilities::PointsKey
1613  p0(nummodes[0],LibUtilities::eGaussGaussLegendre);
1614  LibUtilities::PointsKey
1615  p1(nummodes[1],LibUtilities::eGaussGaussLegendre);
1616  LibUtilities::PointsKey t0(
1617  m_base[0]->GetNumModes(),
1618  LibUtilities::eGaussGaussLegendre);
1619  LibUtilities::PointsKey t1(
1620  m_base[1]->GetNumModes(),
1621  LibUtilities::eGaussGaussLegendre);
1622  LibUtilities::Interp2D(p0, p1, data, t0, t1, coeffs);
1623  }
1624  break;
1625  default:
1626  ASSERTL0(false,
1627  "basis is either not set up or not hierarchicial");
1628  }
1629  }
1630 
1631 
1632  StdRegions::Orientation QuadExp::v_GetTraceOrient(int edge)
1633  {
1634  return m_geom->GetEorient(edge);
1635  }
1636 
1637 
1638  const LibUtilities::BasisSharedPtr& QuadExp::v_GetBasis(int dir) const
1639  {
1640  ASSERTL1(dir >= 0 &&dir <= 1, "input dir is out of range");
1641  return m_base[dir];
1642  }
1643 
1644 
1645  int QuadExp::v_GetNumPoints(const int dir) const
1646  {
1647  return GetNumPoints(dir);
1648  }
1649 
1650  DNekMatSharedPtr QuadExp::v_GenMatrix(
1651  const StdRegions::StdMatrixKey &mkey)
1652  {
1653  DNekMatSharedPtr returnval;
1654  switch (mkey.GetMatrixType())
1655  {
1656  case StdRegions::eHybridDGHelmholtz:
1657  case StdRegions::eHybridDGLamToU:
1658  case StdRegions::eHybridDGLamToQ0:
1659  case StdRegions::eHybridDGLamToQ1:
1660  case StdRegions::eHybridDGLamToQ2:
1661  case StdRegions::eHybridDGHelmBndLam:
1662  case StdRegions::eInvLaplacianWithUnityMean:
1663  returnval = Expansion2D::v_GenMatrix(mkey);
1664  break;
1665  default:
1666  returnval = StdQuadExp::v_GenMatrix(mkey);
1667  }
1668  return returnval;
1669  }
1670 
1671  DNekMatSharedPtr QuadExp::v_CreateStdMatrix(
1672  const StdRegions::StdMatrixKey &mkey)
1673  {
1674  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
1675  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
1676  StdRegions::StdQuadExpSharedPtr tmp =
1677  MemoryManager<StdQuadExp>::AllocateSharedPtr(bkey0,bkey1);
1678  return tmp->GetStdMatrix(mkey);
1679  }
1680 
1681 
1682  DNekScalMatSharedPtr QuadExp::CreateMatrix(const MatrixKey &mkey)
1683  {
1684  DNekScalMatSharedPtr returnval;
1685  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
1686 
1687  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,
1688  "Geometric information is not set up");
1689 
1690  switch (mkey.GetMatrixType())
1691  {
1692  case StdRegions::eMass:
1693  {
1694  if ((m_metricinfo->GetGtype() ==
1695  SpatialDomains::eDeformed) || (mkey.GetNVarCoeff()))
1696  {
1697  NekDouble one = 1.0;
1698  DNekMatSharedPtr mat = GenMatrix(mkey);
1699  returnval = MemoryManager<DNekScalMat>::
1700  AllocateSharedPtr(one,mat);
1701  }
1702  else
1703  {
1704  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1705  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1706  returnval = MemoryManager<DNekScalMat>::
1707  AllocateSharedPtr(jac,mat);
1708  }
1709  }
1710  break;
1711  case StdRegions::eInvMass:
1712  {
1713  if ((m_metricinfo->GetGtype() ==
1714  SpatialDomains::eDeformed) || (mkey.GetNVarCoeff()))
1715  {
1716  NekDouble one = 1.0;
1717  StdRegions::StdMatrixKey masskey(
1718  StdRegions::eMass, DetShapeType(), *this);
1719  DNekMatSharedPtr mat = GenMatrix(masskey);
1720  mat->Invert();
1721 
1722  returnval = MemoryManager<DNekScalMat>::
1723  AllocateSharedPtr(one,mat);
1724  }
1725  else
1726  {
1727  NekDouble fac = 1.0/(m_metricinfo->GetJac(ptsKeys))[0];
1728  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1729  returnval = MemoryManager<DNekScalMat>::
1730  AllocateSharedPtr(fac,mat);
1731  }
1732  }
1733  break;
1734  case StdRegions::eWeakDeriv0:
1735  case StdRegions::eWeakDeriv1:
1736  case StdRegions::eWeakDeriv2:
1737  {
1738  if((m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1739  || (mkey.GetNVarCoeff()))
1740  {
1741  NekDouble one = 1.0;
1742  DNekMatSharedPtr mat = GenMatrix(mkey);
1743 
1744  returnval = MemoryManager<DNekScalMat>::
1745  AllocateSharedPtr(one,mat);
1746  }
1747  else
1748  {
1749  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1750  Array<TwoD, const NekDouble> df =
1751  m_metricinfo->GetDerivFactors(ptsKeys);
1752  int dir = 0;
1753 
1754  switch(mkey.GetMatrixType())
1755  {
1756  case StdRegions::eWeakDeriv0:
1757  dir = 0;
1758  break;
1759  case StdRegions::eWeakDeriv1:
1760  dir = 1;
1761  break;
1762  case StdRegions::eWeakDeriv2:
1763  dir = 2;
1764  break;
1765  default:
1766  break;
1767  }
1768 
1769  MatrixKey deriv0key(StdRegions::eWeakDeriv0,
1770  mkey.GetShapeType(), *this);
1771  MatrixKey deriv1key(StdRegions::eWeakDeriv1,
1772  mkey.GetShapeType(), *this);
1773 
1774  DNekMat &deriv0 = *GetStdMatrix(deriv0key);
1775  DNekMat &deriv1 = *GetStdMatrix(deriv1key);
1776 
1777  int rows = deriv0.GetRows();
1778  int cols = deriv1.GetColumns();
1779 
1780  DNekMatSharedPtr WeakDeriv = MemoryManager<DNekMat>::
1781  AllocateSharedPtr(rows,cols);
1782  (*WeakDeriv) = df[2*dir][0]*deriv0 +
1783  df[2*dir+1][0]*deriv1;
1784  returnval = MemoryManager<DNekScalMat>::
1785  AllocateSharedPtr(jac,WeakDeriv);
1786  }
1787  }
1788  break;
1789  case StdRegions::eLaplacian:
1790  {
1791  if( (m_metricinfo->GetGtype() ==
1792  SpatialDomains::eDeformed) || (mkey.GetNVarCoeff() > 0)
1793  || (mkey.ConstFactorExists
1794  (StdRegions::eFactorSVVCutoffRatio)))
1795  {
1796  NekDouble one = 1.0;
1797  DNekMatSharedPtr mat = GenMatrix(mkey);
1798 
1799  returnval = MemoryManager<DNekScalMat>::
1800  AllocateSharedPtr(one,mat);
1801  }
1802  else
1803  {
1804  MatrixKey lap00key(StdRegions::eLaplacian00,
1805  mkey.GetShapeType(), *this);
1806  MatrixKey lap01key(StdRegions::eLaplacian01,
1807  mkey.GetShapeType(), *this);
1808  MatrixKey lap11key(StdRegions::eLaplacian11,
1809  mkey.GetShapeType(), *this);
1810 
1811  DNekMat &lap00 = *GetStdMatrix(lap00key);
1812  DNekMat &lap01 = *GetStdMatrix(lap01key);
1813  DNekMat &lap11 = *GetStdMatrix(lap11key);
1814 
1815  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1816  Array<TwoD, const NekDouble>
1817  gmat = m_metricinfo->GetGmat(ptsKeys);
1818 
1819  int rows = lap00.GetRows();
1820  int cols = lap00.GetColumns();
1821 
1822  DNekMatSharedPtr lap =
1823  MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols);
1824 
1825  (*lap) = gmat[0][0] * lap00 +
1826  gmat[1][0] * (lap01 + Transpose(lap01)) +
1827  gmat[3][0] * lap11;
1828 
1829  returnval = MemoryManager<DNekScalMat>::
1830  AllocateSharedPtr(jac,lap);
1831  }
1832  }
1833  break;
1834  case StdRegions::eInvLaplacianWithUnityMean:
1835  {
1836  DNekMatSharedPtr mat = GenMatrix(mkey);
1837  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0,mat);
1838  }
1839  break;
1840  case StdRegions::eHelmholtz:
1841  {
1842  NekDouble lambda =
1843  mkey.GetConstFactor(StdRegions::eFactorLambda);
1844 
1845  MatrixKey masskey(mkey, StdRegions::eMass);
1846  DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
1847 
1848  MatrixKey lapkey(mkey, StdRegions::eLaplacian);
1849  DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
1850 
1851  int rows = LapMat.GetRows();
1852  int cols = LapMat.GetColumns();
1853 
1854  DNekMatSharedPtr helm = MemoryManager<DNekMat>::
1855  AllocateSharedPtr(rows,cols);
1856 
1857  NekDouble one = 1.0;
1858  (*helm) = LapMat + lambda*MassMat;
1859 
1860  returnval =
1861  MemoryManager<DNekScalMat>::AllocateSharedPtr(one,helm);
1862  }
1863  break;
1864  case StdRegions::eIProductWRTBase:
1865  {
1866  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1867  {
1868  NekDouble one = 1.0;
1869  DNekMatSharedPtr mat = GenMatrix(mkey);
1870  returnval = MemoryManager<DNekScalMat>::
1871  AllocateSharedPtr(one,mat);
1872  }
1873  else
1874  {
1875  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1876  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1877  returnval = MemoryManager<DNekScalMat>::
1878  AllocateSharedPtr(jac,mat);
1879  }
1880  }
1881  break;
1882  case StdRegions::eIProductWRTDerivBase0:
1883  case StdRegions::eIProductWRTDerivBase1:
1884  case StdRegions::eIProductWRTDerivBase2:
1885  {
1886  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1887  {
1888  NekDouble one = 1.0;
1889  DNekMatSharedPtr mat = GenMatrix(mkey);
1890  returnval = MemoryManager<DNekScalMat>::
1891  AllocateSharedPtr(one,mat);
1892  }
1893  else
1894  {
1895  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1896  const Array<TwoD, const NekDouble>& df =
1897  m_metricinfo->GetDerivFactors(ptsKeys);
1898  int dir = 0;
1899 
1900  switch(mkey.GetMatrixType())
1901  {
1902  case StdRegions::eIProductWRTDerivBase0:
1903  dir = 0;
1904  break;
1905  case StdRegions::eIProductWRTDerivBase1:
1906  dir = 1;
1907  break;
1908  case StdRegions::eIProductWRTDerivBase2:
1909  dir = 2;
1910  break;
1911  default:
1912  break;
1913  }
1914 
1915  MatrixKey iProdDeriv0Key(
1916  StdRegions::eIProductWRTDerivBase0,
1917  mkey.GetShapeType(), *this);
1918  MatrixKey iProdDeriv1Key(
1919  StdRegions::eIProductWRTDerivBase1,
1920  mkey.GetShapeType(), *this);
1921 
1922  DNekMat &stdiprod0 = *GetStdMatrix(iProdDeriv0Key);
1923  DNekMat &stdiprod1 = *GetStdMatrix(iProdDeriv0Key);
1924 
1925  int rows = stdiprod0.GetRows();
1926  int cols = stdiprod1.GetColumns();
1927 
1928  DNekMatSharedPtr mat = MemoryManager<DNekMat>::
1929  AllocateSharedPtr(rows,cols);
1930  (*mat) = df[2*dir][0]*stdiprod0 +
1931  df[2*dir+1][0]*stdiprod1;
1932 
1933  returnval = MemoryManager<DNekScalMat>::
1934  AllocateSharedPtr(jac,mat);
1935  }
1936  }
1937  break;
1938  case StdRegions::eInvHybridDGHelmholtz:
1939  {
1940  NekDouble one = 1.0;
1941 
1942  MatrixKey hkey(StdRegions::eHybridDGHelmholtz,
1943  DetShapeType(), *this,
1944  mkey.GetConstFactors(), mkey.GetVarCoeffs());
1945  DNekMatSharedPtr mat = GenMatrix(hkey);
1946 
1947  mat->Invert();
1948  returnval =
1949  MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1950  }
1951  break;
1952  case StdRegions::eInterpGauss:
1953  {
1954  DNekMatSharedPtr m_Ix;
1955  Array<OneD, NekDouble> coords(1, 0.0);
1956  StdRegions::ConstFactorMap factors = mkey.GetConstFactors();
1957  int edge = (int)factors[StdRegions::eFactorGaussEdge];
1958 
1959  coords[0] = (edge == 0 || edge == 3) ? -1.0 : 1.0;
1960 
1961  m_Ix = m_base[(edge + 1) % 2]->GetI(coords);
1962  returnval =
1963  MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0,m_Ix);
1964  }
1965  break;
1966  case StdRegions::ePreconLinearSpace:
1967  {
1968  NekDouble one = 1.0;
1969  MatrixKey helmkey(
1970  StdRegions::eHelmholtz, mkey.GetShapeType(), *this,
1971  mkey.GetConstFactors(), mkey.GetVarCoeffs());
1972  DNekScalBlkMatSharedPtr helmStatCond =
1973  GetLocStaticCondMatrix(helmkey);
1974  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1975  DNekMatSharedPtr R=BuildVertexMatrix(A);
1976 
1977  returnval =
1978  MemoryManager<DNekScalMat>::AllocateSharedPtr(one, R);
1979  }
1980  break;
1981  default:
1982  {
1983  NekDouble one = 1.0;
1984  DNekMatSharedPtr mat = GenMatrix(mkey);
1985 
1986  returnval =
1987  MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
1988  }
1989  break;
1990  }
1991 
1992  return returnval;
1993  }
1994 
1995 
1996  DNekScalBlkMatSharedPtr QuadExp::CreateStaticCondMatrix(
1997  const MatrixKey &mkey)
1998  {
1999  DNekScalBlkMatSharedPtr returnval;
2000 
2001  ASSERTL2(m_metricinfo->GetGtype()
2002  != SpatialDomains::eNoGeomType,
2003  "Geometric information is not set up");
2004 
2005  // set up block matrix system
2006  unsigned int nbdry = NumBndryCoeffs();
2007  unsigned int nint = (unsigned int)(m_ncoeffs - nbdry);
2008  unsigned int exp_size[] = {nbdry,nint};
2009  unsigned int nblks = 2;
2010  returnval = MemoryManager<DNekScalBlkMat>::
2011  AllocateSharedPtr(nblks,nblks,exp_size,exp_size);
2012  //Really need a constructor which takes Arrays
2013  NekDouble factor = 1.0;
2014 
2015  switch (mkey.GetMatrixType())
2016  {
2017  // this can only use stdregions statically condensed system
2018  // for mass matrix
2019  case StdRegions::eMass:
2020  if ((m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
2021  ||(mkey.GetNVarCoeff()))
2022  {
2023  factor = 1.0;
2024  goto UseLocRegionsMatrix;
2025  }
2026  else
2027  {
2028  factor = (m_metricinfo->GetJac(GetPointsKeys()))[0];
2029  goto UseStdRegionsMatrix;
2030  }
2031  break;
2032  default: // use Deformed case for both
2033  // regular and deformed geometries
2034  factor = 1.0;
2035  goto UseLocRegionsMatrix;
2036  break;
2037  UseStdRegionsMatrix:
2038  {
2039  NekDouble invfactor = 1.0/factor;
2040  NekDouble one = 1.0;
2041  DNekBlkMatSharedPtr mat = GetStdStaticCondMatrix(mkey);
2042  DNekScalMatSharedPtr Atmp;
2043  DNekMatSharedPtr Asubmat;
2044 
2045  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::
2046  AllocateSharedPtr(factor,Asubmat = mat->GetBlock(0,0)));
2047  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::
2048  AllocateSharedPtr(one,Asubmat = mat->GetBlock(0,1)));
2049  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::
2050  AllocateSharedPtr(factor,Asubmat = mat->GetBlock(1,0)));
2051  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::
2052  AllocateSharedPtr(invfactor,Asubmat = mat->GetBlock(1,1)));
2053  }
2054  break;
2055  UseLocRegionsMatrix:
2056  {
2057  int i,j;
2058  NekDouble invfactor = 1.0/factor;
2059  NekDouble one = 1.0;
2060  DNekScalMat &mat = *GetLocMatrix(mkey);
2061  DNekMatSharedPtr A = MemoryManager<DNekMat>::
2062  AllocateSharedPtr(nbdry,nbdry);
2063  DNekMatSharedPtr B = MemoryManager<DNekMat>::
2064  AllocateSharedPtr(nbdry,nint);
2065  DNekMatSharedPtr C = MemoryManager<DNekMat>::
2066  AllocateSharedPtr(nint,nbdry);
2067  DNekMatSharedPtr D = MemoryManager<DNekMat>::
2068  AllocateSharedPtr(nint,nint);
2069 
2070  Array<OneD,unsigned int> bmap(nbdry);
2071  Array<OneD,unsigned int> imap(nint);
2072  GetBoundaryMap(bmap);
2073  GetInteriorMap(imap);
2074 
2075  for (i = 0; i < nbdry; ++i)
2076  {
2077  for(j = 0; j < nbdry; ++j)
2078  {
2079  (*A)(i,j) = mat(bmap[i],bmap[j]);
2080  }
2081 
2082  for(j = 0; j < nint; ++j)
2083  {
2084  (*B)(i,j) = mat(bmap[i],imap[j]);
2085  }
2086  }
2087 
2088  for (i = 0; i < nint; ++i)
2089  {
2090  for(j = 0; j < nbdry; ++j)
2091  {
2092  (*C)(i,j) = mat(imap[i],bmap[j]);
2093  }
2094 
2095  for(j = 0; j < nint; ++j)
2096  {
2097  (*D)(i,j) = mat(imap[i],imap[j]);
2098  }
2099  }
2100 
2101  // Calculate static condensed system
2102  if(nint)
2103  {
2104  D->Invert();
2105  (*B) = (*B)*(*D);
2106  (*A) = (*A) - (*B)*(*C);
2107  }
2108 
2109  DNekScalMatSharedPtr Atmp;
2110 
2111  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::
2112  AllocateSharedPtr(factor, A));
2113  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::
2114  AllocateSharedPtr(one, B));
2115  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::
2116  AllocateSharedPtr(factor, C));
2117  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::
2118  AllocateSharedPtr(invfactor, D));
2119 
2120  }
2121  }
2122  return returnval;
2123  }
2124 
2125 
2126  DNekScalMatSharedPtr QuadExp::v_GetLocMatrix(const MatrixKey &mkey)
2127  {
2128  return m_matrixManager[mkey];
2129  }
2130 
2131 
2132  DNekScalBlkMatSharedPtr QuadExp::v_GetLocStaticCondMatrix(
2133  const MatrixKey &mkey)
2134  {
2135  return m_staticCondMatrixManager[mkey];
2136  }
2137 
2138  void QuadExp::v_DropLocStaticCondMatrix(const MatrixKey &mkey)
2139  {
2140  m_staticCondMatrixManager.DeleteObject(mkey);
2141  }
2142 
2143 
2144  void QuadExp::v_MassMatrixOp(
2145  const Array<OneD, const NekDouble> &inarray,
2146  Array<OneD,NekDouble> &outarray,
2147  const StdRegions::StdMatrixKey &mkey)
2148  {
2149  StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
2150  }
2151 
2152 
2153  void QuadExp::v_LaplacianMatrixOp(
2154  const Array<OneD, const NekDouble> &inarray,
2155  Array<OneD,NekDouble> &outarray,
2156  const StdRegions::StdMatrixKey &mkey)
2157  {
2158  QuadExp::LaplacianMatrixOp_MatFree(inarray, outarray, mkey);
2159  }
2160 
2161 
2162  void QuadExp::v_LaplacianMatrixOp(
2163  const int k1,
2164  const int k2,
2165  const Array<OneD, const NekDouble> &inarray,
2166  Array<OneD,NekDouble> &outarray,
2167  const StdRegions::StdMatrixKey &mkey)
2168  {
2169  StdExpansion::LaplacianMatrixOp_MatFree(
2170  k1, k2, inarray, outarray, mkey);
2171  }
2172 
2173 
2174  void QuadExp::v_WeakDerivMatrixOp(
2175  const int i,
2176  const Array<OneD, const NekDouble> &inarray,
2177  Array<OneD,NekDouble> &outarray,
2178  const StdRegions::StdMatrixKey &mkey)
2179  {
2180  StdExpansion::WeakDerivMatrixOp_MatFree(i, inarray, outarray, mkey);
2181  }
2182 
2183 
2184  void QuadExp::v_WeakDirectionalDerivMatrixOp(
2185  const Array<OneD, const NekDouble> &inarray,
2186  Array<OneD,NekDouble> &outarray,
2187  const StdRegions::StdMatrixKey &mkey)
2188  {
2189  StdExpansion::WeakDirectionalDerivMatrixOp_MatFree(
2190  inarray, outarray, mkey);
2191  }
2192 
2193 
2194  void QuadExp::v_MassLevelCurvatureMatrixOp(
2195  const Array<OneD, const NekDouble> &inarray,
2196  Array<OneD,NekDouble> &outarray,
2197  const StdRegions::StdMatrixKey &mkey)
2198  {
2199  StdExpansion::MassLevelCurvatureMatrixOp_MatFree(
2200  inarray, outarray, mkey);
2201  }
2202 
2203 
2204  void QuadExp::v_HelmholtzMatrixOp(
2205  const Array<OneD, const NekDouble> &inarray,
2206  Array<OneD,NekDouble> &outarray,
2207  const StdRegions::StdMatrixKey &mkey)
2208  {
2209  QuadExp::HelmholtzMatrixOp_MatFree(inarray, outarray, mkey);
2210  }
2211 
2212 
2213  void QuadExp::v_GeneralMatrixOp_MatOp(
2214  const Array<OneD, const NekDouble> &inarray,
2215  Array<OneD,NekDouble> &outarray,
2216  const StdRegions::StdMatrixKey &mkey)
2217  {
2218  MatrixKey newkey(mkey);
2219  DNekScalMatSharedPtr mat = GetLocMatrix(newkey);
2220 
2221  if (inarray.get() == outarray.get())
2222  {
2223  Array<OneD,NekDouble> tmp(m_ncoeffs);
2224  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
2225 
2226  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs, mat->Scale(),
2227  (mat->GetOwnedMatrix())->GetPtr().get(), m_ncoeffs,
2228  tmp.get(), 1, 0.0, outarray.get(), 1);
2229  }
2230  else
2231  {
2232  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),
2233  (mat->GetOwnedMatrix())->GetPtr().get(), m_ncoeffs,
2234  inarray.get(), 1, 0.0, outarray.get(), 1);
2235  }
2236  }
2237 
2238  void QuadExp::v_ReduceOrderCoeffs(
2239  int numMin,
2240  const Array<OneD, const NekDouble> &inarray,
2241  Array<OneD, NekDouble> &outarray)
2242  {
2243  int n_coeffs = inarray.size();
2244 
2245  Array<OneD, NekDouble> coeff (n_coeffs);
2246  Array<OneD, NekDouble> coeff_tmp(n_coeffs, 0.0);
2247  Array<OneD, NekDouble> tmp, tmp2;
2248 
2249  int nmodes0 = m_base[0]->GetNumModes();
2250  int nmodes1 = m_base[1]->GetNumModes();
2251  int numMax = nmodes0;
2252 
2253  Vmath::Vcopy(n_coeffs,inarray,1,coeff_tmp,1);
2254 
2255  const LibUtilities::PointsKey Pkey0(
2256  nmodes0, LibUtilities::eGaussLobattoLegendre);
2257  const LibUtilities::PointsKey Pkey1(
2258  nmodes1, LibUtilities::eGaussLobattoLegendre);
2259  LibUtilities::BasisKey b0(
2260  m_base[0]->GetBasisType(), nmodes0, Pkey0);
2261  LibUtilities::BasisKey b1(
2262  m_base[1]->GetBasisType(), nmodes1, Pkey1);
2263  LibUtilities::BasisKey bortho0(
2264  LibUtilities::eOrtho_A, nmodes0, Pkey0);
2265  LibUtilities::BasisKey bortho1(
2266  LibUtilities::eOrtho_A, nmodes1, Pkey1);
2267 
2268  LibUtilities::InterpCoeff2D(
2269  b0, b1, coeff_tmp, bortho0, bortho1, coeff);
2270 
2271  Vmath::Zero(n_coeffs, coeff_tmp, 1);
2272 
2273  int cnt = 0;
2274  for (int i = 0; i < numMin+1; ++i)
2275  {
2276  Vmath::Vcopy(numMin,
2277  tmp = coeff+cnt,1,
2278  tmp2 = coeff_tmp+cnt,1);
2279 
2280  cnt = i*numMax;
2281  }
2282 
2283  LibUtilities::InterpCoeff2D(
2284  bortho0, bortho1, coeff_tmp,
2285  b0, b1, outarray);
2286  }
2287 
2288  void QuadExp::v_LaplacianMatrixOp_MatFree_Kernel(
2289  const Array<OneD, const NekDouble> &inarray,
2290  Array<OneD, NekDouble> &outarray,
2291  Array<OneD, NekDouble> &wsp)
2292  {
2293  if (m_metrics.count(eMetricLaplacian00) == 0)
2294  {
2295  ComputeLaplacianMetric();
2296  }
2297 
2298  int nquad0 = m_base[0]->GetNumPoints();
2299  int nquad1 = m_base[1]->GetNumPoints();
2300  int nqtot = nquad0*nquad1;
2301  int nmodes0 = m_base[0]->GetNumModes();
2302  int nmodes1 = m_base[1]->GetNumModes();
2303  int wspsize = max(max(max(nqtot,m_ncoeffs),nquad1*nmodes0),nquad0*nmodes1);
2304 
2305  ASSERTL1(wsp.size() >= 3*wspsize,
2306  "Workspace is of insufficient size.");
2307 
2308  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
2309  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
2310  const Array<OneD, const NekDouble>& dbase0 = m_base[0]->GetDbdata();
2311  const Array<OneD, const NekDouble>& dbase1 = m_base[1]->GetDbdata();
2312  const Array<OneD, const NekDouble>& metric00 = m_metrics[eMetricLaplacian00];
2313  const Array<OneD, const NekDouble>& metric01 = m_metrics[eMetricLaplacian01];
2314  const Array<OneD, const NekDouble>& metric11 = m_metrics[eMetricLaplacian11];
2315 
2316  // Allocate temporary storage
2317  Array<OneD,NekDouble> wsp0(wsp);
2318  Array<OneD,NekDouble> wsp1(wsp+wspsize);
2319  Array<OneD,NekDouble> wsp2(wsp+2*wspsize);
2320 
2321  StdExpansion2D::PhysTensorDeriv(inarray,wsp1,wsp2);
2322 
2323  // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
2324  // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
2325  // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
2326  // especially for this purpose
2327  Vmath::Vvtvvtp(nqtot,&metric00[0],1,&wsp1[0],1,&metric01[0],1,&wsp2[0],1,&wsp0[0],1);
2328  Vmath::Vvtvvtp(nqtot,&metric01[0],1,&wsp1[0],1,&metric11[0],1,&wsp2[0],1,&wsp2[0],1);
2329 
2330  // outarray = m = (D_xi1 * B)^T * k
2331  // wsp1 = n = (D_xi2 * B)^T * l
2332  IProductWRTBase_SumFacKernel(dbase0,base1,wsp0,outarray,wsp1,false,true);
2333  IProductWRTBase_SumFacKernel(base0,dbase1,wsp2,wsp1, wsp0,true,false);
2334 
2335  // outarray = outarray + wsp1
2336  // = L * u_hat
2337  Vmath::Vadd(m_ncoeffs,wsp1.get(),1,outarray.get(),1,outarray.get(),1);
2338  }
2339 
2340  void QuadExp::v_ComputeLaplacianMetric()
2341  {
2342  if (m_metrics.count(eMetricQuadrature) == 0)
2343  {
2344  ComputeQuadratureMetric();
2345  }
2346 
2347  const SpatialDomains::GeomType type = m_metricinfo->GetGtype();
2348  const unsigned int nqtot = GetTotPoints();
2349  const unsigned int dim = 2;
2350  const MetricType m[3][3] = { {eMetricLaplacian00, eMetricLaplacian01, eMetricLaplacian02},
2351  {eMetricLaplacian01, eMetricLaplacian11, eMetricLaplacian12},
2352  {eMetricLaplacian02, eMetricLaplacian12, eMetricLaplacian22}
2353  };
2354 
2355  const Array<TwoD, const NekDouble> gmat =
2356  m_metricinfo->GetGmat(GetPointsKeys());
2357  for (unsigned int i = 0; i < dim; ++i)
2358  {
2359  for (unsigned int j = i; j < dim; ++j)
2360  {
2361  m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
2362  if (type == SpatialDomains::eDeformed)
2363  {
2364  Vmath::Vcopy(nqtot, &gmat[i*dim+j][0], 1,
2365  &m_metrics[m[i][j]][0], 1);
2366  }
2367  else
2368  {
2369  Vmath::Fill(nqtot, gmat[i*dim+j][0],
2370  &m_metrics[m[i][j]][0], 1);
2371  }
2372  MultiplyByQuadratureMetric(m_metrics[m[i][j]],
2373  m_metrics[m[i][j]]);
2374 
2375  }
2376  }
2377  }
2378 
2379  void QuadExp::v_SVVLaplacianFilter(
2380  Array<OneD, NekDouble> &array,
2381  const StdRegions::StdMatrixKey &mkey)
2382  {
2383  int nq = GetTotPoints();
2384 
2385  // Calculate sqrt of the Jacobian
2386  Array<OneD, const NekDouble> jac =
2387  m_metricinfo->GetJac(GetPointsKeys());
2388  Array<OneD, NekDouble> sqrt_jac(nq);
2389  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
2390  {
2391  Vmath::Vsqrt(nq,jac,1,sqrt_jac,1);
2392  }
2393  else
2394  {
2395  Vmath::Fill(nq,sqrt(jac[0]),sqrt_jac,1);
2396  }
2397 
2398  // Multiply array by sqrt(Jac)
2399  Vmath::Vmul(nq,sqrt_jac,1,array,1,array,1);
2400 
2401  // Apply std region filter
2402  StdQuadExp::v_SVVLaplacianFilter( array, mkey);
2403 
2404  // Divide by sqrt(Jac)
2405  Vmath::Vdiv(nq,array,1,sqrt_jac,1,array,1);
2406  }
2407 
2408  }//end of namespace
2409 }//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:15
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::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: Expansion2D.cpp:680
Nektar::LocalRegions::Expansion2D::GetGeom2D
SpatialDomains::Geometry2DSharedPtr GetGeom2D() const
Definition: Expansion2D.h:204
Nektar::LocalRegions::Expansion2D::m_edgeNormals
std::map< int, NormalVector > m_edgeNormals
Definition: Expansion2D.h:121
Nektar::LocalRegions::Expansion::m_elmtBndNormDirElmtLen
std::map< int, Array< OneD, NekDouble > > m_elmtBndNormDirElmtLen
the element length in each element boundary(Vertex, edge or face) normal direction calculated based o...
Definition: Expansion.h:284
Nektar::LocalRegions::Expansion::GetGeom
SpatialDomains::GeometrySharedPtr GetGeom() const
Definition: Expansion.cpp:172
Nektar::LocalRegions::Expansion::m_geom
SpatialDomains::GeometrySharedPtr m_geom
Definition: Expansion.h:272
Nektar::LocalRegions::Expansion::GetLeftAdjacentElementExp
ExpansionSharedPtr GetLeftAdjacentElementExp() const
Definition: Expansion.h:438
Nektar::LocalRegions::Expansion::BuildVertexMatrix
DNekMatSharedPtr BuildVertexMatrix(const DNekScalMatSharedPtr &r_bnd)
Definition: Expansion.cpp:103
Nektar::LocalRegions::Expansion::GetLeftAdjacentElementTrace
int GetLeftAdjacentElementTrace() const
Definition: Expansion.h:454
Nektar::LocalRegions::Expansion::m_metricinfo
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
Definition: Expansion.h:273
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:318
Nektar::LocalRegions::Expansion::GetTraceOrient
StdRegions::Orientation GetTraceOrient(int trace)
Definition: Expansion.h:164
Nektar::LocalRegions::Expansion::GetLocMatrix
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition: Expansion.cpp:90
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:705
Nektar::LocalRegions::QuadExp::v_GetLocMatrix
virtual DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey)
Definition: QuadExp.cpp:2126
Nektar::LocalRegions::QuadExp::v_GetMetricInfo
virtual const SpatialDomains::GeomFactorsSharedPtr & v_GetMetricInfo() const
Definition: QuadExp.cpp:1521
Nektar::LocalRegions::QuadExp::CreateStaticCondMatrix
DNekScalBlkMatSharedPtr CreateStaticCondMatrix(const MatrixKey &mkey)
Definition: QuadExp.cpp:1996
Nektar::LocalRegions::QuadExp::v_ComputeTraceNormal
virtual void v_ComputeTraceNormal(const int edge)
Definition: QuadExp.cpp:1220
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:569
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:2379
Nektar::LocalRegions::QuadExp::v_CreateStdMatrix
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: QuadExp.cpp:1671
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:2174
Nektar::LocalRegions::QuadExp::v_ReduceOrderCoeffs
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: QuadExp.cpp:2238
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:2144
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:2153
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:697
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:2204
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:1650
Nektar::LocalRegions::QuadExp::v_GetTraceOrient
virtual StdRegions::Orientation v_GetTraceOrient(int edge)
Definition: QuadExp.cpp:1632
Nektar::LocalRegions::QuadExp::v_GetTracePhysMap
virtual void v_GetTracePhysMap(const int edge, Array< OneD, int > &outarray)
Definition: QuadExp.cpp:920
Nektar::LocalRegions::QuadExp::v_AlignVectorToCollapsedDir
virtual void v_AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: QuadExp.cpp:524
Nektar::LocalRegions::QuadExp::CreateMatrix
DNekScalMatSharedPtr CreateMatrix(const MatrixKey &mkey)
Definition: QuadExp.cpp:1682
Nektar::LocalRegions::QuadExp::v_GetEdgePhysVals
virtual void v_GetEdgePhysVals(const int edge, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: QuadExp.cpp:721
Nektar::LocalRegions::QuadExp::v_GetCoord
virtual void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords)
Definition: QuadExp.cpp:674
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:2340
Nektar::LocalRegions::QuadExp::v_GetLocStaticCondMatrix
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey)
Definition: QuadExp.cpp:2132
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:665
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:610
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:786
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:2194
Nektar::LocalRegions::QuadExp::v_DropLocStaticCondMatrix
void v_DropLocStaticCondMatrix(const MatrixKey &mkey)
Definition: QuadExp.cpp:2138
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:2213
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:2184
Nektar::LocalRegions::QuadExp::v_GetTraceQFactors
virtual void v_GetTraceQFactors(const int edge, Array< OneD, NekDouble > &outarray)
Definition: QuadExp.cpp:968
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:2288
Nektar::LocalRegions::QuadExp::v_GetBasis
virtual const LibUtilities::BasisSharedPtr & v_GetBasis(int dir) const
Definition: QuadExp.cpp:1638
Nektar::LocalRegions::QuadExp::v_GetStdExp
virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const
Definition: QuadExp.cpp:646
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:1533
Nektar::LocalRegions::QuadExp::v_GetLinStdExp
virtual StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const
Definition: QuadExp.cpp:654
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:1645
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:854
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:171
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:210
Nektar::StdRegions::StdExpansion::GetBoundaryMap
void GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:687
Nektar::StdRegions::StdExpansion::GetNcoeffs
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:124
Nektar::StdRegions::StdExpansion::GetTotPoints
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:134
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:158
Nektar::StdRegions::StdExpansion::HelmholtzMatrixOp_MatFree
void HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:1288
Nektar::StdRegions::StdExpansion::GetStdMatrix
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:617
Nektar::StdRegions::StdExpansion::GetPointsKeys
const LibUtilities::PointsKeyVector GetPointsKeys() const
Definition: StdExpansion.h:1078
Nektar::StdRegions::StdExpansion::GetStdStaticCondMatrix
DNekBlkMatSharedPtr GetStdStaticCondMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:622
Nektar::StdRegions::StdExpansion::MassMatrixOp
void MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:762
Nektar::StdRegions::StdExpansion::GetLocStaticCondMatrix
DNekScalBlkMatSharedPtr GetLocStaticCondMatrix(const LocalRegions::MatrixKey &mkey)
Definition: StdExpansion.h:660
Nektar::StdRegions::StdExpansion::IProductWRTDerivBase_SumFac
void IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:1203
Nektar::StdRegions::StdExpansion::NormVectorIProductWRTBase
void NormVectorIProductWRTBase(const Array< OneD, const NekDouble > &Fx, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:628
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:537
Nektar::StdRegions::StdExpansion::LaplacianMatrixOp_MatFree
void LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:1241
Nektar::StdRegions::StdExpansion::GetTraceToElementMap
void GetTraceToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation traceOrient=eForwards, int P=-1, int Q=-1)
Definition: StdExpansion.h:703
Nektar::StdRegions::StdExpansion::DetShapeType
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:376
Nektar::StdRegions::StdExpansion::GetInteriorMap
void GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:692
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:433
Nektar::StdRegions::StdExpansion::GenMatrix
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:850
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:208
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:1745
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:221
Nektar::StdRegions::StdExpansion::IProductWRTBase_SumFac
void IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Definition: StdExpansion.h:1147
Nektar::StdRegions::StdExpansion::MultiplyByQuadratureMetric
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:733
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition: StdExpansion.h:1162
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:265
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:197
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:167
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:272
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:273
Nektar::StdRegions::StdExpansionSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1714
Nektar::StdRegions::ConstFactorMap
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:314
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
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:820
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:475
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:691
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:192
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:565
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:513
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:322
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:225
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:291
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:257
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:436
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:1226
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:625
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1199
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:372
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:267
Generated on Wed Nov 24 2021 10:33:28 for Nektar++ by doxygen 1.9.1