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