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