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