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