Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
TriExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File TriExp.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 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Expasion for triangular elements.
32 //
33 ///////////////////////////////////////////////////////////////////////////////
35 #include <LocalRegions/TriExp.h>
36 #include <LocalRegions/SegExp.h>
40 
41 using namespace std;
42 
43 namespace Nektar
44 {
45  namespace LocalRegions
46  {
47  TriExp::TriExp(const LibUtilities::BasisKey &Ba,
48  const LibUtilities::BasisKey &Bb,
50  StdExpansion (LibUtilities::StdTriData::getNumberOfCoefficients(Ba.GetNumModes(),(Bb.GetNumModes())),2,Ba,Bb),
51  StdExpansion2D(LibUtilities::StdTriData::getNumberOfCoefficients(Ba.GetNumModes(),(Bb.GetNumModes())),Ba,Bb),
52  StdTriExp(Ba,Bb),
53  Expansion (geom),
54  Expansion2D (geom),
55  m_matrixManager(
56  boost::bind(&TriExp::CreateMatrix, this, _1),
57  std::string("TriExpMatrix")),
58  m_staticCondMatrixManager(
59  boost::bind(&TriExp::CreateStaticCondMatrix, this, _1),
60  std::string("TriExpStaticCondMatrix"))
61  {
62  }
63 
64 
66  StdExpansion(T),
67  StdExpansion2D(T),
68  StdTriExp(T),
69  Expansion(T),
70  Expansion2D(T),
71  m_matrixManager(T.m_matrixManager),
72  m_staticCondMatrixManager(T.m_staticCondMatrixManager)
73  {
74  }
75 
76 
78  {
79  }
80 
81 
82 
84  {
85  int nquad0 = m_base[0]->GetNumPoints();
86  int nquad1 = m_base[1]->GetNumPoints();
88  NekDouble ival;
89  Array<OneD,NekDouble> tmp(nquad0*nquad1);
90 
91  // multiply inarray with Jacobian
92  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
93  {
94  Vmath::Vmul(nquad0*nquad1, jac, 1, inarray, 1,tmp, 1);
95  }
96  else
97  {
98  Vmath::Smul(nquad0*nquad1, jac[0], inarray, 1, tmp, 1);
99  }
100 
101  // call StdQuadExp version;
102  ival = StdTriExp::v_Integral(tmp);
103  return ival;
104  }
105 
106 
108  Array<OneD,NekDouble> &out_d0,
109  Array<OneD,NekDouble> &out_d1,
110  Array<OneD,NekDouble> &out_d2)
111  {
112  int nquad0 = m_base[0]->GetNumPoints();
113  int nquad1 = m_base[1]->GetNumPoints();
114  int nqtot = nquad0*nquad1;
116  = m_metricinfo->GetDerivFactors(GetPointsKeys());
117 
118  Array<OneD,NekDouble> diff0(2*nqtot);
119  Array<OneD,NekDouble> diff1(diff0+nqtot);
120 
121  StdTriExp::v_PhysDeriv(inarray, diff0, diff1);
122 
123  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
124  {
125  if(out_d0.num_elements())
126  {
127  Vmath::Vmul (nqtot,df[0],1,diff0,1, out_d0, 1);
128  Vmath::Vvtvp (nqtot,df[1],1,diff1,1, out_d0, 1, out_d0,1);
129  }
130 
131  if(out_d1.num_elements())
132  {
133  Vmath::Vmul (nqtot,df[2],1,diff0,1, out_d1, 1);
134  Vmath::Vvtvp (nqtot,df[3],1,diff1,1, out_d1, 1, out_d1,1);
135  }
136 
137  if(out_d2.num_elements())
138  {
139  Vmath::Vmul (nqtot,df[4],1,diff0,1, out_d2, 1);
140  Vmath::Vvtvp (nqtot,df[5],1,diff1,1, out_d2, 1, out_d2,1);
141  }
142  }
143  else // regular geometry
144  {
145  if(out_d0.num_elements())
146  {
147  Vmath::Smul (nqtot, df[0][0], diff0, 1, out_d0, 1);
148  Blas::Daxpy (nqtot, df[1][0], diff1, 1, out_d0, 1);
149  }
150 
151  if(out_d1.num_elements())
152  {
153  Vmath::Smul (nqtot, df[2][0], diff0, 1, out_d1, 1);
154  Blas::Daxpy (nqtot, df[3][0], diff1, 1, out_d1, 1);
155  }
156 
157  if(out_d2.num_elements())
158  {
159  Vmath::Smul (nqtot, df[4][0], diff0, 1, out_d2, 1);
160  Blas::Daxpy (nqtot, df[5][0], diff1, 1, out_d2, 1);
161  }
162  }
163  }
164 
165 
166  void TriExp::v_PhysDeriv(const int dir,
167  const Array<OneD, const NekDouble>& inarray,
168  Array<OneD, NekDouble> &outarray)
169  {
170  switch(dir)
171  {
172  case 0:
173  {
175  }
176  break;
177  case 1:
178  {
180  }
181  break;
182  case 2:
183  {
185  }
186  break;
187  default:
188  {
189  ASSERTL1(false,"input dir is out of range");
190  }
191  break;
192  }
193  }
194 
196  const Array<OneD, const NekDouble> &inarray,
197  const Array<OneD, const NekDouble> &direction,
199  {
200  if(! out.num_elements())
201  {
202  return;
203  }
204 
205  int nquad0 = m_base[0]->GetNumPoints();
206  int nquad1 = m_base[1]->GetNumPoints();
207  int nqtot = nquad0*nquad1;
208 
209  const Array<TwoD, const NekDouble>& df =
210  m_metricinfo->GetDerivFactors(GetPointsKeys());
211 
212  Array<OneD,NekDouble> diff0(2*nqtot);
213  Array<OneD,NekDouble> diff1(diff0+nqtot);
214 
215  // diff0 = du/d_xi, diff1 = du/d_eta
216  StdTriExp::v_PhysDeriv(inarray, diff0, diff1);
217 
218  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
219  {
221 
222 
223  // D^v_xi = v_x*d_xi/dx + v_y*d_xi/dy + v_z*d_xi/dz
224  // D^v_eta = v_x*d_eta/dx + v_y*d_eta/dy + v_z*d_eta/dz
225  for (int i=0; i< 2; ++i)
226  {
227  tangmat[i] = Array<OneD, NekDouble>(nqtot,0.0);
228  for (int k=0; k<(m_geom->GetCoordim()); ++k)
229  {
230  Vmath::Vvtvp(nqtot,&df[2*k+i][0],1,&direction[k*nqtot],1,&tangmat[i][0],1,&tangmat[i][0],1);
231  }
232  }
233 
234  /// D_v = D^v_xi * du/d_xi + D^v_eta * du/d_eta
235  Vmath::Vmul (nqtot,&tangmat[0][0],1,&diff0[0],1, &out[0], 1);
236  Vmath::Vvtvp (nqtot,&tangmat[1][0],1,&diff1[0],1, &out[0], 1, &out[0],1);
237  }
238  else
239  {
240  ASSERTL1(m_metricinfo->GetGtype() == SpatialDomains::eDeformed,"Wrong route");
241  }
242  }
243 
244 
246  Array<OneD,NekDouble> &outarray)
247  {
248  IProductWRTBase(inarray,outarray);
249 
250  // get Mass matrix inverse
252  DetShapeType(),*this);
253  DNekScalMatSharedPtr matsys = m_matrixManager[masskey];
254 
255  // copy inarray in case inarray == outarray
256  NekVector<NekDouble> in (m_ncoeffs,outarray,eCopy);
258 
259  out = (*matsys)*in;
260  }
261 
262 
264  Array<OneD, NekDouble> &outarray)
265  {
266  int i,j;
267  int npoints[2] = {m_base[0]->GetNumPoints(),
268  m_base[1]->GetNumPoints()};
269  int nmodes[2] = {m_base[0]->GetNumModes(),
270  m_base[1]->GetNumModes()};
271 
272  fill(outarray.get(), outarray.get()+m_ncoeffs, 0.0 );
273 
274  Array<OneD, NekDouble> physEdge[3];
275  Array<OneD, NekDouble> coeffEdge[3];
276  for(i = 0; i < 3; i++)
277  {
278  // define physEdge and add 1 so can interpolate grl10 points if necessary
279  physEdge[i] = Array<OneD, NekDouble>(max(npoints[i!=0],npoints[0]));
280  coeffEdge[i] = Array<OneD, NekDouble>(nmodes[i!=0]);
281  }
282 
283  for(i = 0; i < npoints[0]; i++)
284  {
285  physEdge[0][i] = inarray[i];
286  }
287 
288  // extract data in cartesian directions
289  for(i = 0; i < npoints[1]; i++)
290  {
291  physEdge[1][i] = inarray[npoints[0]-1+i*npoints[0]];
292  physEdge[2][i] = inarray[i*npoints[0]];
293  }
294 
295  SegExpSharedPtr segexp[3];
296  segexp[0] = MemoryManager<LocalRegions::SegExp>::AllocateSharedPtr(m_base[0]->GetBasisKey(),GetGeom2D()->GetEdge(0));
297 
299  {
300  for(i = 1; i < 3; i++)
301  {
302  segexp[i] = MemoryManager<LocalRegions::SegExp>::AllocateSharedPtr(m_base[i!=0]->GetBasisKey(),GetGeom2D()->GetEdge(i));
303  }
304  }
305  else // interploate using edge 0 GLL distribution
306  {
307  for(i = 1; i < 3; i++)
308  {
309  segexp[i] = MemoryManager<LocalRegions::SegExp>::AllocateSharedPtr(m_base[0]->GetBasisKey(),GetGeom2D()->GetEdge(i));
310 
311  LibUtilities::Interp1D(m_base[1]->GetPointsKey(),physEdge[i],
312  m_base[0]->GetPointsKey(),physEdge[i]);
313  }
314  npoints[1] = npoints[0];
315  }
316 
317 
318  Array<OneD, unsigned int> mapArray;
319  Array<OneD, int> signArray;
320  NekDouble sign;
321  // define an orientation to get EdgeToElmtMapping from Cartesian data
324 
325  for(i = 0; i < 3; i++)
326  {
327  segexp[i]->FwdTrans_BndConstrained(physEdge[i],coeffEdge[i]);
328 
329  // this orient goes with the one above and so could
330  // probably set both to eForwards
331  GetEdgeToElementMap(i,orient[i],mapArray,signArray);
332  for(j=0; j < nmodes[i!=0]; j++)
333  {
334  sign = (NekDouble) signArray[j];
335  outarray[ mapArray[j] ] = sign * coeffEdge[i][j];
336  }
337  }
338 
339  int nBoundaryDofs = NumBndryCoeffs();
340  int nInteriorDofs = m_ncoeffs - nBoundaryDofs;
341 
342  if (nInteriorDofs > 0) {
345 
347  MassMatrixOp(outarray,tmp0,stdmasskey);
348  IProductWRTBase(inarray,tmp1);
349 
350  Vmath::Vsub(m_ncoeffs, tmp1, 1, tmp0, 1, tmp1, 1);
351 
352  // get Mass matrix inverse (only of interior DOF)
353  // use block (1,1) of the static condensed system
354  // note: this block alreay contains the inverse matrix
355  MatrixKey masskey(StdRegions::eMass,DetShapeType(),*this);
356  DNekScalMatSharedPtr matsys = (m_staticCondMatrixManager[masskey])->GetBlock(1,1);
357 
358  Array<OneD, NekDouble> rhs(nInteriorDofs);
359  Array<OneD, NekDouble> result(nInteriorDofs);
360 
361  GetInteriorMap(mapArray);
362 
363  for(i = 0; i < nInteriorDofs; i++)
364  {
365  rhs[i] = tmp1[ mapArray[i] ];
366  }
367 
368  Blas::Dgemv('N', nInteriorDofs, nInteriorDofs, matsys->Scale(), &((matsys->GetOwnedMatrix())->GetPtr())[0],
369  nInteriorDofs,rhs.get(),1,0.0,result.get(),1);
370 
371  for(i = 0; i < nInteriorDofs; i++)
372  {
373  outarray[ mapArray[i] ] = result[i];
374  }
375  }
376  }
377 
378 
380  Array<OneD, NekDouble> &outarray)
381  {
382  IProductWRTBase_SumFac(inarray,outarray);
383  }
384 
385 
386  void TriExp::v_IProductWRTDerivBase(const int dir,
387  const Array<OneD, const NekDouble>& inarray,
388  Array<OneD, NekDouble> & outarray)
389  {
390  IProductWRTDerivBase_SumFac(dir,inarray,outarray);
391  }
392 
393 
395  Array<OneD, NekDouble> &outarray,
396  bool multiplybyweights)
397  {
398  int nquad0 = m_base[0]->GetNumPoints();
399  int nquad1 = m_base[1]->GetNumPoints();
400  int order0 = m_base[0]->GetNumModes();
401 
402  if(multiplybyweights)
403  {
404  Array<OneD,NekDouble> tmp(nquad0*nquad1+nquad1*order0);
405  Array<OneD,NekDouble> wsp(tmp+nquad0*nquad1);
406 
407  MultiplyByQuadratureMetric(inarray,tmp);
408  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),m_base[1]->GetBdata(),tmp,outarray,wsp);
409  }
410  else
411  {
412  Array<OneD,NekDouble> wsp(+nquad1*order0);
413 
414  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),m_base[1]->GetBdata(),
415  inarray,outarray,wsp);
416  }
417  }
418 
419 
421  Array<OneD, NekDouble> &outarray)
422  {
423  int nq = GetTotPoints();
425  DNekScalMatSharedPtr iprodmat = m_matrixManager[iprodmatkey];
426 
427  Blas::Dgemv('N',m_ncoeffs,nq,iprodmat->Scale(),(iprodmat->GetOwnedMatrix())->GetPtr().get(),
428  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
429 
430  }
431 
432 
434  const Array<OneD, const NekDouble>& inarray,
435  Array<OneD, NekDouble> & outarray)
436  {
437  ASSERTL1((dir==0)||(dir==1)||(dir==2),"Invalid direction.");
438  ASSERTL1((dir==2)?(m_geom->GetCoordim()==3):true,"Invalid direction.");
439 
440  int i;
441  int nquad0 = m_base[0]->GetNumPoints();
442  int nquad1 = m_base[1]->GetNumPoints();
443  int nqtot = nquad0*nquad1;
444  int nmodes0 = m_base[0]->GetNumModes();
445  int wspsize = max(max(nqtot,m_ncoeffs),nquad1*nmodes0);
446 
447  const Array<TwoD, const NekDouble>& df =
448  m_metricinfo->GetDerivFactors(GetPointsKeys());
449 
450  Array<OneD, NekDouble> tmp0 (6*wspsize);
451  Array<OneD, NekDouble> tmp1 (tmp0 + wspsize);
452  Array<OneD, NekDouble> tmp2 (tmp0 + 2*wspsize);
453  Array<OneD, NekDouble> tmp3 (tmp0 + 3*wspsize);
454  Array<OneD, NekDouble> gfac0(tmp0 + 4*wspsize);
455  Array<OneD, NekDouble> gfac1(tmp0 + 5*wspsize);
456 
457  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
458  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
459 
460  // set up geometric factor: 2/(1-z1)
461  for(i = 0; i < nquad1; ++i)
462  {
463  gfac0[i] = 2.0/(1-z1[i]);
464  }
465  for(i = 0; i < nquad0; ++i)
466  {
467  gfac1[i] = 0.5*(1+z0[i]);
468  }
469 
470  for(i = 0; i < nquad1; ++i)
471  {
472  Vmath::Smul(nquad0,gfac0[i],&inarray[0]+i*nquad0,1,&tmp0[0]+i*nquad0,1);
473  }
474 
475  for(i = 0; i < nquad1; ++i)
476  {
477  Vmath::Vmul(nquad0,&gfac1[0],1,&tmp0[0]+i*nquad0,1,&tmp1[0]+i*nquad0,1);
478  }
479 
480  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
481  {
482  Vmath::Vmul(nqtot,&df[2*dir][0], 1,&tmp0[0], 1,&tmp0[0],1);
483  Vmath::Vmul(nqtot,&df[2*dir+1][0],1,&tmp1[0], 1,&tmp1[0],1);
484  Vmath::Vmul(nqtot,&df[2*dir+1][0],1,&inarray[0],1,&tmp2[0],1);
485  }
486  else
487  {
488  Vmath::Smul(nqtot, df[2*dir][0], tmp0, 1, tmp0, 1);
489  Vmath::Smul(nqtot, df[2*dir+1][0], tmp1, 1, tmp1, 1);
490  Vmath::Smul(nqtot, df[2*dir+1][0], inarray, 1, tmp2, 1);
491  }
492  Vmath::Vadd(nqtot, tmp0, 1, tmp1, 1, tmp1, 1);
493 
494  MultiplyByQuadratureMetric(tmp1,tmp1);
495  MultiplyByQuadratureMetric(tmp2,tmp2);
496 
497  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),m_base[1]->GetBdata() ,tmp1,tmp3 ,tmp0);
498  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata() ,m_base[1]->GetDbdata(),tmp2,outarray,tmp0);
499  Vmath::Vadd(m_ncoeffs, tmp3, 1, outarray, 1, outarray, 1);
500  }
501 
502 
504  const Array<OneD, const NekDouble>& inarray,
505  Array<OneD, NekDouble> &outarray)
506  {
507  int nq = GetTotPoints();
509 
510  switch(dir)
511  {
512  case 0:
513  {
515  }
516  break;
517  case 1:
518  {
520  }
521  break;
522  case 2:
523  {
525  }
526  break;
527  default:
528  {
529  ASSERTL1(false,"input dir is out of range");
530  }
531  break;
532  }
533 
534  MatrixKey iprodmatkey(mtype,DetShapeType(),*this);
535  DNekScalMatSharedPtr iprodmat = m_matrixManager[iprodmatkey];
536 
537  Blas::Dgemv('N',m_ncoeffs,nq,iprodmat->Scale(),(iprodmat->GetOwnedMatrix())->GetPtr().get(),
538  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
539 
540  }
541 
546  Array<OneD, NekDouble> &outarray)
547  {
548  int nq = m_base[0]->GetNumPoints()*m_base[1]->GetNumPoints();
550 
551  const Array<OneD, const Array<OneD, NekDouble> > &normals =
552  GetLeftAdjacentElementExp()->GetFaceNormal(
554 
555  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
556  {
557  Vmath::Vvtvvtp(nq,&normals[0][0],1,&Fx[0],1,
558  &normals[1][0],1,&Fy[0],1,&Fn[0],1);
559  Vmath::Vvtvp (nq,&normals[2][0],1,&Fz[0],1,&Fn[0],1,&Fn[0],1);
560  }
561  else
562  {
563  Vmath::Svtsvtp(nq,normals[0][0],&Fx[0],1,
564  normals[1][0],&Fy[0],1,&Fn[0],1);
565  Vmath::Svtvp (nq,normals[2][0],&Fz[0],1,&Fn[0],1,&Fn[0],1);
566  }
567 
568  IProductWRTBase(Fn,outarray);
569  }
570 
572  const Array<OneD, const Array<OneD, NekDouble> > &Fvec,
573  Array<OneD, NekDouble> &outarray)
574  {
575  NormVectorIProductWRTBase(Fvec[0], Fvec[1], Fvec[2], outarray);
576  }
577 
579  {
580 
582  ::AllocateSharedPtr(m_base[0]->GetBasisKey(),
583  m_base[1]->GetBasisKey());
584  }
585 
587  Array<OneD,NekDouble> &coords)
588  {
589  int i;
590 
591  ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[1] <= 1.0 &&
592  Lcoords[1] >= -1.0 && Lcoords[1] <=1.0,
593  "Local coordinates are not in region [-1,1]");
594 
595  m_geom->FillGeom();
596 
597  for(i = 0; i < m_geom->GetCoordim(); ++i)
598  {
599  coords[i] = m_geom->GetCoord(i,Lcoords);
600  }
601  }
602 
604  Array<OneD, NekDouble> &coords_0,
605  Array<OneD, NekDouble> &coords_1,
606  Array<OneD, NekDouble> &coords_2)
607  {
608  Expansion::v_GetCoords(coords_0, coords_1, coords_2);
609  }
610 
611 
612  /**
613  * Given the local cartesian coordinate \a Lcoord evaluate the
614  * value of physvals at this point by calling through to the
615  * StdExpansion method
616  */
618  const Array<OneD, const NekDouble> &Lcoord,
619  const Array<OneD, const NekDouble> &physvals)
620  {
621  // Evaluate point in local (eta) coordinates.
622  return StdTriExp::v_PhysEvaluate(Lcoord,physvals);
623  }
624 
626  {
628 
629  ASSERTL0(m_geom,"m_geom not defined");
630  m_geom->GetLocCoords(coord,Lcoord);
631 
632  return StdTriExp::v_PhysEvaluate(Lcoord, physvals);
633  }
634 
635 
637  const int edge,
638  const StdRegions::StdExpansionSharedPtr &EdgeExp,
639  const Array<OneD, const NekDouble> &inarray,
640  Array<OneD,NekDouble> &outarray,
642  {
643  v_GetEdgePhysVals(edge,EdgeExp,inarray,outarray);
644  }
645 
647  const int edge,
648  const Array<OneD, const NekDouble> &inarray,
649  Array<OneD,NekDouble> &outarray)
650  {
651  int nquad0 = m_base[0]->GetNumPoints();
652  int nquad1 = m_base[1]->GetNumPoints();
653 
654  StdRegions::Orientation edgedir = GetEorient(edge);
655  switch(edge)
656  {
657  case 0:
658  if (edgedir == StdRegions::eForwards)
659  {
660  Vmath::Vcopy(nquad0,&(inarray[0]),1,&(outarray[0]),1);
661  }
662  else
663  {
664  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0-1),-1,
665  &(outarray[0]),1);
666  }
667  break;
668  case 1:
669  if (edgedir == StdRegions::eForwards)
670  {
671  Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0-1),nquad0,
672  &(outarray[0]),1);
673  }
674  else
675  {
676  Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0*nquad1-1),
677  -nquad0, &(outarray[0]),1);
678  }
679  break;
680  case 2:
681  if (edgedir == StdRegions::eForwards)
682  {
683  Vmath::Vcopy(nquad1,&(inarray[0]) + nquad0*(nquad1-1),
684  -nquad0,&(outarray[0]),1);
685  }
686  else
687  {
688  Vmath::Vcopy(nquad1,&(inarray[0]),nquad0,
689  &(outarray[0]),1);
690  }
691  break;
692  default:
693  ASSERTL0(false,"edge value (< 3) is out of range");
694  break;
695  }
696  }
697 
699  const Array<OneD, const NekDouble> &inarray,
700  Array<OneD,NekDouble> &outarray)
701  {
702  int nquad0 = m_base[0]->GetNumPoints();
703  int nquad1 = m_base[1]->GetNumPoints();
704 
705  // get points in Cartesian orientation
706  switch(edge)
707  {
708  case 0:
709  Vmath::Vcopy(nquad0, &(inarray[0]), 1, &(outarray[0]), 1);
710  break;
711  case 1:
712  Vmath::Vcopy(nquad1, &(inarray[0])+(nquad0-1),
713  nquad0, &(outarray[0]), 1);
714  break;
715  case 2:
716  Vmath::Vcopy(nquad1, &(inarray[0]), nquad0, &(outarray[0]), 1);
717  break;
718  default:
719  ASSERTL0(false,"edge value (< 3) is out of range");
720  break;
721  }
722 
723  // Interpolate if required
724  if(m_base[edge?1:0]->GetPointsKey() != EdgeExp->GetBasis(0)->GetPointsKey())
725  {
726  Array<OneD,NekDouble> outtmp(max(nquad0,nquad1));
727 
728  outtmp = outarray;
729 
730  LibUtilities::Interp1D(m_base[edge?1:0]->GetPointsKey(),
731  outtmp,
732  EdgeExp->GetBasis(0)->GetPointsKey(),
733  outarray);
734  }
735 
736  //Reverse data if necessary
738  {
739  Vmath::Reverse(EdgeExp->GetNumPoints(0),&outarray[0],1,
740  &outarray[0],1);
741  }
742 
743  }
744 
745 
747  const int edge,const Array<OneD, const NekDouble> &inarray,
748  Array<OneD, NekDouble> &outarray)
749  {
750  ASSERTL0(false,
751  "Routine not implemented for triangular elements");
752  }
753 
755  const int edge,
756  Array<OneD, NekDouble> &outarray)
757  {
758  ASSERTL0(false,
759  "Routine not implemented for triangular elements");
760  }
761 
762 
763 
765  const int edge,
766  Array<OneD, int> &outarray)
767  {
768  int nquad0 = m_base[0]->GetNumPoints();
769  int nquad1 = m_base[1]->GetNumPoints();
770 
771  // Get points in Cartesian orientation
772  switch (edge)
773  {
774  case 0:
775  outarray = Array<OneD, int>(nquad0);
776  for (int i = 0; i < nquad0; ++i)
777  {
778  outarray[i] = i;
779  }
780  break;
781  case 1:
782  outarray = Array<OneD, int>(nquad1);
783  for (int i = 0; i < nquad1; ++i)
784  {
785  outarray[i] = (nquad0-1) + i * nquad0;
786  }
787  break;
788  case 2:
789  outarray = Array<OneD, int>(nquad1);
790  for (int i = 0; i < nquad1; ++i)
791  {
792  outarray[i] = i*nquad0;
793  }
794  break;
795  default:
796  ASSERTL0(false, "edge value (< 3) is out of range");
797  break;
798  }
799 
800  }
801 
802 
803  void TriExp::v_ComputeEdgeNormal(const int edge)
804  {
805  int i;
806  const SpatialDomains::GeomFactorsSharedPtr & geomFactors = GetGeom()->GetMetricInfo();
808  const SpatialDomains::GeomType type = geomFactors->GetGtype();
809  const Array<TwoD, const NekDouble> & df = geomFactors->GetDerivFactors(ptsKeys);
810  const Array<OneD, const NekDouble> & jac = geomFactors->GetJac(ptsKeys);
811  int nqe = m_base[0]->GetNumPoints();
812  int dim = GetCoordim();
813 
816  for (i = 0; i < dim; ++i)
817  {
818  normal[i] = Array<OneD, NekDouble>(nqe);
819  }
820 
821  // Regular geometry case
823  {
824  NekDouble fac;
825  // Set up normals
826  switch(edge)
827  {
828  case 0:
829  for(i = 0; i < GetCoordim(); ++i)
830  {
831  Vmath::Fill(nqe,-df[2*i+1][0],normal[i],1);
832  }
833  break;
834  case 1:
835  for(i = 0; i < GetCoordim(); ++i)
836  {
837  Vmath::Fill(nqe,df[2*i+1][0] + df[2*i][0],normal[i],1);
838  }
839  break;
840  case 2:
841  for(i = 0; i < GetCoordim(); ++i)
842  {
843  Vmath::Fill(nqe,-df[2*i][0],normal[i],1);
844  }
845  break;
846  default:
847  ASSERTL0(false,"Edge is out of range (edge < 3)");
848  }
849 
850  // normalise
851  fac = 0.0;
852  for(i =0 ; i < GetCoordim(); ++i)
853  {
854  fac += normal[i][0]*normal[i][0];
855  }
856  fac = 1.0/sqrt(fac);
857  for (i = 0; i < GetCoordim(); ++i)
858  {
859  Vmath::Smul(nqe,fac,normal[i],1,normal[i],1);
860  }
861  }
862  else // Set up deformed normals
863  {
864  int j;
865 
866  int nquad0 = ptsKeys[0].GetNumPoints();
867  int nquad1 = ptsKeys[1].GetNumPoints();
868 
869  LibUtilities::PointsKey from_key;
870 
871  Array<OneD,NekDouble> normals(GetCoordim()*max(nquad0,nquad1),0.0);
872  Array<OneD,NekDouble> edgejac(GetCoordim()*max(nquad0,nquad1),0.0);
873 
874  // Extract Jacobian along edges and recover local
875  // derivates (dx/dr) for polynomial interpolation by
876  // multiplying m_gmat by jacobian
877  switch(edge)
878  {
879  case 0:
880  for(j = 0; j < nquad0; ++j)
881  {
882  edgejac[j] = jac[j];
883  for(i = 0; i < GetCoordim(); ++i)
884  {
885  normals[i*nquad0+j] = -df[2*i+1][j]*edgejac[j];
886  }
887  }
888  from_key = ptsKeys[0];
889  break;
890  case 1:
891  for(j = 0; j < nquad1; ++j)
892  {
893  edgejac[j] = jac[nquad0*j+nquad0-1];
894  for(i = 0; i < GetCoordim(); ++i)
895  {
896  normals[i*nquad1+j] = (df[2*i][nquad0*j + nquad0-1] + df[2*i+1][nquad0*j + nquad0-1])*edgejac[j];
897  }
898  }
899  from_key = ptsKeys[1];
900  break;
901  case 2:
902  for(j = 0; j < nquad1; ++j)
903  {
904  edgejac[j] = jac[nquad0*j];
905  for(i = 0; i < GetCoordim(); ++i)
906  {
907  normals[i*nquad1+j] = -df[2*i][nquad0*j]*edgejac[j];
908  }
909  }
910  from_key = ptsKeys[1];
911  break;
912  default:
913  ASSERTL0(false,"edge is out of range (edge < 3)");
914 
915  }
916 
917  int nq = from_key.GetNumPoints();
918  Array<OneD,NekDouble> work(nqe,0.0);
919 
920  // interpolate Jacobian and invert
921  LibUtilities::Interp1D(from_key,jac,m_base[0]->GetPointsKey(),work);
922  Vmath::Sdiv(nq,1.0,&work[0],1,&work[0],1);
923 
924  // interpolate
925  for(i = 0; i < GetCoordim(); ++i)
926  {
927  LibUtilities::Interp1D(from_key,&normals[i*nq],m_base[0]->GetPointsKey(),&normal[i][0]);
928  Vmath::Vmul(nqe,work,1,normal[i],1,normal[i],1);
929  }
930 
931  //normalise normal vectors
932  Vmath::Zero(nqe,work,1);
933  for(i = 0; i < GetCoordim(); ++i)
934  {
935  Vmath::Vvtvp(nqe,normal[i],1, normal[i],1,work,1,work,1);
936  }
937 
938  Vmath::Vsqrt(nqe,work,1,work,1);
939  Vmath::Sdiv(nqe,1.0,work,1,work,1);
940 
941  for(i = 0; i < GetCoordim(); ++i)
942  {
943  Vmath::Vmul(nqe,normal[i],1,work,1,normal[i],1);
944  }
945 
946  // Reverse direction so that points are in
947  // anticlockwise direction if edge >=2
948  if(edge >= 2)
949  {
950  for(i = 0; i < GetCoordim(); ++i)
951  {
952  Vmath::Reverse(nqe,normal[i],1, normal[i],1);
953  }
954  }
955  }
957  {
958  for(i = 0; i < GetCoordim(); ++i)
959  {
960  if(geomFactors->GetGtype() == SpatialDomains::eDeformed)
961  {
962  Vmath::Reverse(nqe, normal[i], 1, normal[i],1);
963  }
964  }
965  }
966  }
967 
969  {
970  return m_geom->GetCoordim();
971  }
972 
973 
975  const std::vector<unsigned int > &nummodes, const int mode_offset, NekDouble * coeffs)
976  {
977  int data_order0 = nummodes[mode_offset];
978  int fillorder0 = min(m_base[0]->GetNumModes(),data_order0);
979  int data_order1 = nummodes[mode_offset+1];
980  int order1 = m_base[1]->GetNumModes();
981  int fillorder1 = min(order1,data_order1);
982 
983  switch(m_base[0]->GetBasisType())
984  {
986  {
987  int i;
988  int cnt = 0;
989  int cnt1 = 0;
990 
992  "Extraction routine not set up for this basis");
993 
994  Vmath::Zero(m_ncoeffs,coeffs,1);
995  for(i = 0; i < fillorder0; ++i)
996  {
997  Vmath::Vcopy(fillorder1-i,&data[cnt],1,&coeffs[cnt1],1);
998  cnt += data_order1-i;
999  cnt1 += order1-i;
1000  }
1001  }
1002  break;
1003  default:
1004  ASSERTL0(false,"basis is either not set up or not hierarchicial");
1005  }
1006  }
1007 
1008 
1010  {
1011  return GetGeom2D()->GetEorient(edge);
1012  }
1013 
1014 
1016  {
1017  return GetGeom2D()->GetCartesianEorient(edge);
1018  }
1019 
1020 
1022  {
1023  ASSERTL1(dir >= 0 &&dir <= 1,"input dir is out of range");
1024  return m_base[dir];
1025  }
1026 
1027 
1028  int TriExp::v_GetNumPoints(const int dir) const
1029  {
1030  return GetNumPoints(dir);
1031  }
1032 
1033 
1035  {
1036  DNekMatSharedPtr returnval;
1037  switch(mkey.GetMatrixType())
1038  {
1046  returnval = Expansion2D::v_GenMatrix(mkey);
1047  break;
1048  default:
1049  returnval = StdTriExp::v_GenMatrix(mkey);
1050  break;
1051  }
1052 
1053  return returnval;
1054  }
1055 
1056 
1058  {
1059  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
1060  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
1062  AllocateSharedPtr(bkey0,bkey1);
1063 
1064  return tmp->GetStdMatrix(mkey);
1065  }
1066 
1067 
1069  {
1070  DNekScalMatSharedPtr returnval;
1072 
1073  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
1074 
1075  switch(mkey.GetMatrixType())
1076  {
1077  case StdRegions::eMass:
1078  {
1079  if((m_metricinfo->GetGtype() == SpatialDomains::eDeformed)||
1080  (mkey.GetNVarCoeff()))
1081  {
1082  NekDouble one = 1.0;
1083  DNekMatSharedPtr mat = GenMatrix(mkey);
1084  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1085  }
1086  else
1087  {
1088  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1089  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1090  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
1091  }
1092  }
1093  break;
1094  case StdRegions::eInvMass:
1095  {
1096  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1097  {
1098  NekDouble one = 1.0;
1100  *this);
1101  DNekMatSharedPtr mat = GenMatrix(masskey);
1102  mat->Invert();
1103 
1104  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1105  }
1106  else
1107  {
1108  NekDouble fac = 1.0/(m_metricinfo->GetJac(ptsKeys))[0];
1109  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1110  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(fac,mat);
1111 
1112  }
1113  }
1114  break;
1118  {
1119  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed || mkey.GetNVarCoeff())
1120  {
1121  NekDouble one = 1.0;
1122  DNekMatSharedPtr mat = GenMatrix(mkey);
1123 
1124  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1125  }
1126  else
1127  {
1128  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1129  Array<TwoD, const NekDouble> df = m_metricinfo->GetDerivFactors(ptsKeys);
1130  int dir = 0;
1131  switch(mkey.GetMatrixType())
1132  {
1134  dir = 0;
1135  break;
1137  dir = 1;
1138  break;
1140  dir = 2;
1141  break;
1142  default:
1143  break;
1144  }
1145 
1147  mkey.GetShapeType(), *this);
1149  mkey.GetShapeType(), *this);
1150 
1151  DNekMat &deriv0 = *GetStdMatrix(deriv0key);
1152  DNekMat &deriv1 = *GetStdMatrix(deriv1key);
1153 
1154  int rows = deriv0.GetRows();
1155  int cols = deriv1.GetColumns();
1156 
1158  (*WeakDeriv) = df[2*dir][0]*deriv0 + df[2*dir+1][0]*deriv1;
1159 
1160  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,WeakDeriv);
1161  }
1162  }
1163  break;
1165  {
1166  if( (m_metricinfo->GetGtype() == SpatialDomains::eDeformed) ||
1168  {
1169  NekDouble one = 1.0;
1170  DNekMatSharedPtr mat = GenMatrix(mkey);
1171 
1172  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1173  }
1174  else
1175  {
1177  mkey.GetShapeType(), *this);
1179  mkey.GetShapeType(), *this);
1181  mkey.GetShapeType(), *this);
1182 
1183  DNekMat &lap00 = *GetStdMatrix(lap00key);
1184  DNekMat &lap01 = *GetStdMatrix(lap01key);
1185  DNekMat &lap11 = *GetStdMatrix(lap11key);
1186 
1187  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1189  m_metricinfo->GetGmat(ptsKeys);
1190 
1191  int rows = lap00.GetRows();
1192  int cols = lap00.GetColumns();
1193 
1195 
1196  (*lap) = gmat[0][0] * lap00 +
1197  gmat[1][0] * (lap01 + Transpose(lap01)) +
1198  gmat[3][0] * lap11;
1199 
1200  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,lap);
1201  }
1202  }
1203  break;
1205  {
1206  DNekMatSharedPtr mat = GenMatrix(mkey);
1207  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0,mat);
1208  }
1209  break;
1211  {
1213 
1214  MatrixKey masskey(mkey, StdRegions::eMass);
1215  DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
1216 
1217  MatrixKey lapkey(mkey, StdRegions::eLaplacian);
1218  DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
1219 
1220  int rows = LapMat.GetRows();
1221  int cols = LapMat.GetColumns();
1222 
1224 
1225  NekDouble one = 1.0;
1226  (*helm) = LapMat + factor*MassMat;
1227 
1228  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,helm);
1229  }
1230  break;
1232  {
1233  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1234  {
1235  NekDouble one = 1.0;
1236  DNekMatSharedPtr mat = GenMatrix(mkey);
1237  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1238  }
1239  else
1240  {
1241  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1242  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1243  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
1244  }
1245  }
1246  break;
1250  {
1251  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1252  {
1253  NekDouble one = 1.0;
1254  DNekMatSharedPtr mat = GenMatrix(mkey);
1255  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1256  }
1257  else
1258  {
1259  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1260 
1261  const Array<TwoD, const NekDouble>& df = m_metricinfo->GetDerivFactors(ptsKeys);
1262  int dir = 0;
1263 
1264  switch(mkey.GetMatrixType())
1265  {
1267  dir = 0;
1268  break;
1270  dir = 1;
1271  break;
1273  dir = 2;
1274  break;
1275  default:
1276  break;
1277  }
1278 
1280  mkey.GetShapeType(), *this);
1282  mkey.GetShapeType(), *this);
1283 
1284  DNekMat &stdiprod0 = *GetStdMatrix(iProdDeriv0Key);
1285  DNekMat &stdiprod1 = *GetStdMatrix(iProdDeriv0Key);
1286 
1287  int rows = stdiprod0.GetRows();
1288  int cols = stdiprod1.GetColumns();
1289 
1291  (*mat) = df[2*dir][0]*stdiprod0 + df[2*dir+1][0]*stdiprod1;
1292 
1293  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
1294  }
1295  }
1296  break;
1297 
1299  {
1300  NekDouble one = 1.0;
1301 
1303 
1304  DNekMatSharedPtr mat = GenMatrix(hkey);
1305 
1306  mat->Invert();
1307  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1308  }
1309  break;
1311  {
1312  NekDouble one = 1.0;
1313  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1314  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1315  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1317 
1319  }
1320  break;
1321  default:
1322  {
1323  NekDouble one = 1.0;
1324  DNekMatSharedPtr mat = GenMatrix(mkey);
1325 
1326  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1327  }
1328  break;
1329  }
1330 
1331  return returnval;
1332  }
1333 
1334 
1336  {
1337  DNekScalBlkMatSharedPtr returnval;
1339 
1340  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
1341 
1342  // set up block matrix system
1343  unsigned int nbdry = NumBndryCoeffs();
1344  unsigned int nint = (unsigned int)(m_ncoeffs - nbdry);
1345  unsigned int exp_size[] = {nbdry,nint};
1346  unsigned int nblks = 2;
1347  returnval = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nblks,nblks,exp_size,exp_size); //Really need a constructor which takes Arrays
1348  NekDouble factor = 1.0;
1349 
1350  switch(mkey.GetMatrixType())
1351  {
1352  // this can only use stdregions statically condensed system for mass matrix
1353  case StdRegions::eMass:
1354  if((m_metricinfo->GetGtype() == SpatialDomains::eDeformed)||(mkey.GetNVarCoeff()))
1355  {
1356  factor = 1.0;
1357  goto UseLocRegionsMatrix;
1358  }
1359  else
1360  {
1361  factor = (m_metricinfo->GetJac(ptsKeys))[0];
1362  goto UseStdRegionsMatrix;
1363  }
1364  break;
1365  default: // use Deformed case for both regular and deformed geometries
1366  factor = 1.0;
1367  goto UseLocRegionsMatrix;
1368  break;
1369  UseStdRegionsMatrix:
1370  {
1371  NekDouble invfactor = 1.0/factor;
1372  NekDouble one = 1.0;
1374  DNekScalMatSharedPtr Atmp;
1375  DNekMatSharedPtr Asubmat;
1376 
1377  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(0,0)));
1378  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,Asubmat = mat->GetBlock(0,1)));
1379  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(1,0)));
1380  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,Asubmat = mat->GetBlock(1,1)));
1381  }
1382  break;
1383 
1384  UseLocRegionsMatrix:
1385  {
1386  int i,j;
1387  NekDouble invfactor = 1.0/factor;
1388  NekDouble one = 1.0;
1389 
1390  DNekScalMat &mat = *GetLocMatrix(mkey);
1391 
1396 
1397  Array<OneD,unsigned int> bmap(nbdry);
1398  Array<OneD,unsigned int> imap(nint);
1399  GetBoundaryMap(bmap);
1400  GetInteriorMap(imap);
1401 
1402  for(i = 0; i < nbdry; ++i)
1403  {
1404  for(j = 0; j < nbdry; ++j)
1405  {
1406  (*A)(i,j) = mat(bmap[i],bmap[j]);
1407  }
1408 
1409  for(j = 0; j < nint; ++j)
1410  {
1411  (*B)(i,j) = mat(bmap[i],imap[j]);
1412  }
1413  }
1414 
1415  for(i = 0; i < nint; ++i)
1416  {
1417  for(j = 0; j < nbdry; ++j)
1418  {
1419  (*C)(i,j) = mat(imap[i],bmap[j]);
1420  }
1421 
1422  for(j = 0; j < nint; ++j)
1423  {
1424  (*D)(i,j) = mat(imap[i],imap[j]);
1425  }
1426  }
1427 
1428  // Calculate static condensed system
1429  if(nint)
1430  {
1431  D->Invert();
1432  (*B) = (*B)*(*D);
1433  (*A) = (*A) - (*B)*(*C);
1434  }
1435 
1436  DNekScalMatSharedPtr Atmp;
1437 
1438  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,A));
1439  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,B));
1440  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,C));
1441  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,D));
1442 
1443  }
1444  }
1445 
1446  return returnval;
1447  }
1448 
1449 
1451  {
1452  return m_matrixManager[mkey];
1453  }
1454 
1455 
1457  {
1458  return m_staticCondMatrixManager[mkey];
1459  }
1460 
1462  {
1463  m_staticCondMatrixManager.DeleteObject(mkey);
1464  }
1465 
1466 
1467 
1469  Array<OneD,NekDouble> &outarray,
1470  const StdRegions::StdMatrixKey &mkey)
1471  {
1472  StdExpansion::MassMatrixOp_MatFree(inarray,outarray,mkey);
1473  }
1474 
1475 
1477  Array<OneD,NekDouble> &outarray,
1478  const StdRegions::StdMatrixKey &mkey)
1479  {
1480  TriExp::LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
1481  }
1482 
1483 
1484  void TriExp::v_LaplacianMatrixOp(const int k1, const int k2,
1485  const Array<OneD, const NekDouble> &inarray,
1486  Array<OneD,NekDouble> &outarray,
1487  const StdRegions::StdMatrixKey &mkey)
1488  {
1489  StdExpansion::LaplacianMatrixOp_MatFree(k1,k2,inarray,outarray,mkey);
1490  }
1491 
1492 
1494  const Array<OneD, const NekDouble> &inarray,
1495  Array<OneD,NekDouble> &outarray,
1496  const StdRegions::StdMatrixKey &mkey)
1497  {
1498  StdExpansion::WeakDerivMatrixOp_MatFree(i,inarray,outarray,mkey);
1499  }
1500 
1501 
1503  Array<OneD,NekDouble> &outarray,
1504  const StdRegions::StdMatrixKey &mkey)
1505  {
1506  StdExpansion::WeakDirectionalDerivMatrixOp_MatFree(inarray,outarray,mkey);
1507  }
1508 
1509 
1511  Array<OneD,NekDouble> &outarray,
1512  const StdRegions::StdMatrixKey &mkey)
1513  {
1514  StdExpansion::MassLevelCurvatureMatrixOp_MatFree(inarray,outarray,mkey);
1515  }
1516 
1517 
1519  Array<OneD,NekDouble> &outarray,
1520  const StdRegions::StdMatrixKey &mkey)
1521  {
1522  TriExp::HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
1523  }
1524 
1525 
1527  Array<OneD,NekDouble> &outarray,
1528  const StdRegions::StdMatrixKey &mkey)
1529  {
1530  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1531 
1532  if(inarray.get() == outarray.get())
1533  {
1535  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
1536 
1537  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1538  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
1539  }
1540  else
1541  {
1542  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1543  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
1544  }
1545  }
1546 
1547 
1549  const Array<OneD, const NekDouble> &inarray,
1550  Array<OneD, NekDouble> &outarray,
1552  {
1553  if (m_metrics.count(eMetricLaplacian00) == 0)
1554  {
1556  }
1557 
1558  int nquad0 = m_base[0]->GetNumPoints();
1559  int nquad1 = m_base[1]->GetNumPoints();
1560  int nqtot = nquad0*nquad1;
1561  int nmodes0 = m_base[0]->GetNumModes();
1562  int nmodes1 = m_base[1]->GetNumModes();
1563  int wspsize = max(max(max(nqtot,m_ncoeffs),nquad1*nmodes0),nquad0*nmodes1);
1564 
1565  ASSERTL1(wsp.num_elements() >= 3*wspsize,
1566  "Workspace is of insufficient size.");
1567 
1568  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
1569  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
1570  const Array<OneD, const NekDouble>& dbase0 = m_base[0]->GetDbdata();
1571  const Array<OneD, const NekDouble>& dbase1 = m_base[1]->GetDbdata();
1575 
1576  // Allocate temporary storage
1577  Array<OneD,NekDouble> wsp0(wsp);
1578  Array<OneD,NekDouble> wsp1(wsp+wspsize);
1579  Array<OneD,NekDouble> wsp2(wsp+2*wspsize);
1580 
1581  StdExpansion2D::PhysTensorDeriv(inarray,wsp1,wsp2);
1582 
1583  // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1584  // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1585  // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
1586  // especially for this purpose
1587  Vmath::Vvtvvtp(nqtot,&metric00[0],1,&wsp1[0],1,&metric01[0],1,&wsp2[0],1,&wsp0[0],1);
1588  Vmath::Vvtvvtp(nqtot,&metric01[0],1,&wsp1[0],1,&metric11[0],1,&wsp2[0],1,&wsp2[0],1);
1589 
1590  // outarray = m = (D_xi1 * B)^T * k
1591  // wsp1 = n = (D_xi2 * B)^T * l
1592  IProductWRTBase_SumFacKernel(dbase0,base1,wsp0,outarray,wsp1);
1593  IProductWRTBase_SumFacKernel(base0,dbase1,wsp2,wsp1, wsp0);
1594 
1595  // outarray = outarray + wsp1
1596  // = L * u_hat
1597  Vmath::Vadd(m_ncoeffs,wsp1.get(),1,outarray.get(),1,outarray.get(),1);
1598  }
1599 
1600 
1602  {
1603  if (m_metrics.count(eMetricQuadrature) == 0)
1604  {
1606  }
1607 
1608  unsigned int i, j;
1609  const SpatialDomains::GeomType type = m_metricinfo->GetGtype();
1610  const unsigned int nqtot = GetTotPoints();
1611  const unsigned int dim = 2;
1615  };
1616 
1617  Array<OneD, NekDouble> dEta_dXi[2] = {Array<OneD, NekDouble>(nqtot,1.0),
1618  Array<OneD, NekDouble>(nqtot,1.0)};
1619 
1620  for (i = 0; i < dim; ++i)
1621  {
1622  for (j = i; j < dim; ++j)
1623  {
1624  m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
1625  }
1626  }
1627 
1628  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
1629  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
1630  const unsigned int nquad0 = m_base[0]->GetNumPoints();
1631  const unsigned int nquad1 = m_base[1]->GetNumPoints();
1632  const Array<TwoD, const NekDouble>& df =
1633  m_metricinfo->GetDerivFactors(GetPointsKeys());
1634 
1635  for(i = 0; i < nquad1; i++)
1636  {
1637  Blas::Dscal(nquad0,2.0/(1-z1[i]),&dEta_dXi[0][0]+i*nquad0,1);
1638  Blas::Dscal(nquad0,2.0/(1-z1[i]),&dEta_dXi[1][0]+i*nquad0,1);
1639  }
1640  for(i = 0; i < nquad0; i++)
1641  {
1642  Blas::Dscal(nquad1,0.5*(1+z0[i]),&dEta_dXi[1][0]+i,nquad0);
1643  }
1644 
1645  Array<OneD, NekDouble> tmp(nqtot);
1646  if((type == SpatialDomains::eRegular ||
1648  {
1649  Vmath::Smul (nqtot,df[0][0],&dEta_dXi[0][0],1,&tmp[0],1);
1650  Vmath::Svtvp(nqtot,df[1][0],&dEta_dXi[1][0],1,&tmp[0],1,&tmp[0],1);
1651 
1652  Vmath::Vmul (nqtot,&tmp[0],1, &tmp[0],1,&m_metrics[eMetricLaplacian00][0],1);
1653  Vmath::Smul (nqtot,df[1][0],&tmp[0],1,&m_metrics[eMetricLaplacian01][0],1);
1654 
1655 
1656  Vmath::Smul (nqtot,df[2][0],&dEta_dXi[0][0],1,&tmp[0],1);
1657  Vmath::Svtvp(nqtot,df[3][0],&dEta_dXi[1][0],1,&tmp[0],1,&tmp[0],1);
1658 
1659  Vmath::Vvtvp(nqtot,&tmp[0],1, &tmp[0],1,&m_metrics[eMetricLaplacian00][0],1,&m_metrics[eMetricLaplacian00][0],1);
1660  Vmath::Svtvp(nqtot,df[3][0],&tmp[0],1,&m_metrics[eMetricLaplacian01][0],1,&m_metrics[eMetricLaplacian01][0],1);
1661 
1662  if(GetCoordim() == 3)
1663  {
1664  Vmath::Smul (nqtot,df[4][0],&dEta_dXi[0][0],1,&tmp[0],1);
1665  Vmath::Svtvp(nqtot,df[5][0],&dEta_dXi[1][0],1,&tmp[0],1,&tmp[0],1);
1666 
1667  Vmath::Vvtvp(nqtot,&tmp[0],1, &tmp[0],1,&m_metrics[eMetricLaplacian00][0],1,&m_metrics[eMetricLaplacian00][0],1);
1668  Vmath::Svtvp(nqtot,df[5][0],&tmp[0],1,&m_metrics[eMetricLaplacian01][0],1,&m_metrics[eMetricLaplacian01][0],1);
1669  }
1670 
1671  NekDouble g2 = df[1][0]*df[1][0] + df[3][0]*df[3][0];
1672  if(GetCoordim() == 3)
1673  {
1674  g2 += df[5][0]*df[5][0];
1675  }
1676  Vmath::Fill(nqtot,g2,&m_metrics[eMetricLaplacian11][0],1);
1677  }
1678  else
1679  {
1680 
1681  Vmath::Vmul (nqtot,&df[0][0],1,&dEta_dXi[0][0],1,&tmp[0],1);
1682  Vmath::Vvtvp(nqtot,&df[1][0],1,&dEta_dXi[1][0],1,&tmp[0],1,&tmp[0],1);
1683 
1684  Vmath::Vmul (nqtot,&tmp[0], 1,&tmp[0], 1,&m_metrics[eMetricLaplacian00][0],1);
1685  Vmath::Vmul (nqtot,&df[1][0],1,&tmp[0], 1,&m_metrics[eMetricLaplacian01][0],1);
1686  Vmath::Vmul (nqtot,&df[1][0],1,&df[1][0],1,&m_metrics[eMetricLaplacian11][0],1);
1687 
1688 
1689  Vmath::Vmul (nqtot,&df[2][0],1,&dEta_dXi[0][0],1,&tmp[0],1);
1690  Vmath::Vvtvp(nqtot,&df[3][0],1,&dEta_dXi[1][0],1,&tmp[0],1,&tmp[0],1);
1691 
1692  Vmath::Vvtvp(nqtot,&tmp[0], 1,&tmp[0], 1,&m_metrics[eMetricLaplacian00][0],1,&m_metrics[eMetricLaplacian00][0],1);
1693  Vmath::Vvtvp(nqtot,&df[3][0],1,&tmp[0], 1,&m_metrics[eMetricLaplacian01][0],1,&m_metrics[eMetricLaplacian01][0],1);
1694  Vmath::Vvtvp(nqtot,&df[3][0],1,&df[3][0],1,&m_metrics[eMetricLaplacian11][0],1,&m_metrics[eMetricLaplacian11][0],1);
1695 
1696  if(GetCoordim() == 3)
1697  {
1698  Vmath::Vmul (nqtot,&df[4][0],1,&dEta_dXi[0][0],1,&tmp[0],1);
1699  Vmath::Vvtvp(nqtot,&df[5][0],1,&dEta_dXi[1][0],1,&tmp[0],1,&tmp[0],1);
1700 
1701  Vmath::Vvtvp(nqtot,&tmp[0], 1,&tmp[0], 1,&m_metrics[eMetricLaplacian00][0],1,&m_metrics[eMetricLaplacian00][0],1);
1702  Vmath::Vvtvp(nqtot,&df[5][0],1,&tmp[0], 1,&m_metrics[eMetricLaplacian01][0],1,&m_metrics[eMetricLaplacian01][0],1);
1703  Vmath::Vvtvp(nqtot,&df[5][0],1,&df[5][0],1,&m_metrics[eMetricLaplacian11][0],1,&m_metrics[eMetricLaplacian11][0],1);
1704  }
1705  }
1706 
1707  for (unsigned int i = 0; i < dim; ++i)
1708  {
1709  for (unsigned int j = i; j < dim; ++j)
1710  {
1712  m_metrics[m[i][j]]);
1713 
1714  }
1715  }
1716  }
1717 
1718  /**
1719  * Function is used to compute exactly the advective numerical flux on
1720  * theinterface of two elements with different expansions, hence an
1721  * appropriate number of Gauss points has to be used. The number of
1722  * Gauss points has to be equal to the number used by the highest
1723  * polynomial degree of the two adjacent elements. Furthermore, this
1724  * function is used to compute the sensor value in each element.
1725  *
1726  * @param numMin Is the reduced polynomial order
1727  * @param inarray Input array of coefficients
1728  * @param dumpVar Output array of reduced coefficients.
1729  */
1731  int numMin,
1732  const Array<OneD, const NekDouble> &inarray,
1733  Array<OneD, NekDouble> &outarray)
1734  {
1735  int n_coeffs = inarray.num_elements();
1736  int nquad0 = m_base[0]->GetNumPoints();
1737  int nquad1 = m_base[1]->GetNumPoints();
1738  int nqtot = nquad0*nquad1;
1739  int nmodes0 = m_base[0]->GetNumModes();
1740  int nmodes1 = m_base[1]->GetNumModes();
1741  int numMin2 = nmodes0, i;
1742 
1743  Array<OneD, NekDouble> coeff(n_coeffs,0.0);
1744  Array<OneD, NekDouble> phys_tmp(nqtot,0.0);
1745  Array<OneD, NekDouble> tmp, tmp2;
1746 
1747  const LibUtilities::PointsKey Pkey0 = m_base[0]->GetPointsKey();
1748  const LibUtilities::PointsKey Pkey1 = m_base[1]->GetPointsKey();
1749 
1751  m_base[0]->GetBasisType(), nmodes0, Pkey0);
1753  m_base[1]->GetBasisType(), nmodes1, Pkey1);
1754  LibUtilities::BasisKey bortho0(
1755  LibUtilities::eOrtho_A, nmodes0, Pkey0);
1756  LibUtilities::BasisKey bortho1(
1757  LibUtilities::eOrtho_B, nmodes1, Pkey1);
1758 
1759  // Check if it is also possible to use the same InterCoeff routine
1760  // which is also used for Quadrilateral and Hexagonal shaped
1761  // elements
1762 
1763  // For now, set up the used basis on the standard element to
1764  // calculate the phys values, set up the orthogonal basis to do a
1765  // forward transform, to obtain the coefficients in orthogonal
1766  // coefficient space
1767  StdRegions::StdTriExpSharedPtr m_OrthoTriExp;
1769 
1771  ::AllocateSharedPtr(b0, b1);
1772  m_OrthoTriExp = MemoryManager<StdRegions::StdTriExp>
1773  ::AllocateSharedPtr(bortho0, bortho1);
1774 
1775  m_TriExp ->BwdTrans(inarray,phys_tmp);
1776  m_OrthoTriExp->FwdTrans(phys_tmp, coeff);
1777 
1778  for (i = 0; i < n_coeffs; i++)
1779  {
1780  if (i == numMin)
1781  {
1782  coeff[i] = 0.0;
1783  numMin += numMin2 - 1;
1784  numMin2 -= 1.0;
1785  }
1786  }
1787 
1788  m_OrthoTriExp->BwdTrans(coeff,phys_tmp);
1789  m_TriExp ->FwdTrans(phys_tmp, outarray);
1790  }
1791 
1793  Array<OneD, NekDouble> &array,
1794  const StdRegions::StdMatrixKey &mkey)
1795  {
1796  int nq = GetTotPoints();
1797 
1798  // Calculate sqrt of the Jacobian
1800  m_metricinfo->GetJac(GetPointsKeys());
1801  Array<OneD, NekDouble> sqrt_jac(nq);
1802  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1803  {
1804  Vmath::Vsqrt(nq,jac,1,sqrt_jac,1);
1805  }
1806  else
1807  {
1808  Vmath::Fill(nq,sqrt(jac[0]),sqrt_jac,1);
1809  }
1810 
1811  // Multiply array by sqrt(Jac)
1812  Vmath::Vmul(nq,sqrt_jac,1,array,1,array,1);
1813 
1814  // Apply std region filter
1815  StdTriExp::v_SVVLaplacianFilter( array, mkey);
1816 
1817  // Divide by sqrt(Jac)
1818  Vmath::Vdiv(nq,array,1,sqrt_jac,1,array,1);
1819  }
1820 
1821  }
1822 }
1823 
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
Definition: TriExp.h:284
const LibUtilities::PointsKeyVector GetPointsKeys() const
boost::shared_ptr< StdTriExp > StdTriExpSharedPtr
Definition: StdTriExp.h:267
virtual NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals)
Definition: TriExp.cpp:617
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:470
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:122
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
StdRegions::Orientation GetCartesianEorient(int edge)
Definition: StdExpansion.h:772
const ConstFactorMap & GetConstFactors() const
Definition: StdMatrixKey.h:142
const VarCoeffMap & GetVarCoeffs() const
Definition: StdMatrixKey.h:168
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:220
virtual int v_GetCoordim()
Definition: TriExp.cpp:968
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:82
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey)
Definition: TriExp.cpp:1456
void MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:971
static Array< OneD, NekDouble > NullNekDouble1DArray
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
Definition: TriExp.h:285
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:22
void IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:394
virtual int v_GetNumPoints(const int dir) const
Definition: TriExp.cpp:1028
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:942
virtual void v_MassLevelCurvatureMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1510
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:386
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
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
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
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: TriExp.cpp:107
void LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
virtual void v_WeakDirectionalDerivMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1502
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
Definition: StdExpansion.h:229
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
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
Principle Modified Functions .
Definition: BasisType.h:49
void NormVectorIProductWRTBase(const Array< OneD, const NekDouble > &Fx, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:727
std::map< int, StdRegions::NormalVector > m_edgeNormals
Definition: Expansion2D.h:135
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
Definition: Expansion.h:126
StdRegions::Orientation GetEorient(int edge)
Definition: StdExpansion.h:762
STL namespace.
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: TriExp.cpp:974
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
DNekMatSharedPtr BuildVertexMatrix(const DNekScalMatSharedPtr &r_bnd)
Definition: Expansion.cpp:98
LibUtilities::ShapeType GetShapeType() const
Definition: StdMatrixKey.h:87
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
Definition: TriExp.cpp:603
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
SpatialDomains::GeometrySharedPtr m_geom
Definition: Expansion.h:125
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: TriExp.cpp:195
DNekScalBlkMatSharedPtr GetLocStaticCondMatrix(const LocalRegions::MatrixKey &mkey)
Definition: StdExpansion.h:747
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1034
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: TriExp.cpp:245
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:700
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
void IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
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: TriExp.cpp:646
virtual void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1518
Principle Orthogonal Functions .
Definition: BasisType.h:47
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:131
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:141
virtual void v_IProductWRTDerivBase_MatOp(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:503
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:433
virtual void v_ComputeEdgeNormal(const int edge)
Definition: TriExp.cpp:803
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: TriExp.cpp:542
DNekBlkMatSharedPtr GetStdStaticCondMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:705
boost::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:266
virtual void v_ComputeLaplacianMetric()
Definition: TriExp.cpp:1601
SpatialDomains::Geometry2DSharedPtr GetGeom2D() const
Definition: Expansion2D.h:269
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1071
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
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[p]*base1[pq] and put into ou...
Definition: TriExp.cpp:379
virtual void v_IProductWRTBase_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:420
virtual NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
Integrates the specified function over the domain.
Definition: TriExp.cpp:83
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: TriExp.cpp:636
Principle Modified Functions .
Definition: BasisType.h:50
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
Definition: Expansion.cpp:213
virtual void v_GetEdgeQFactors(const int edge, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:754
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
virtual DNekScalBlkMatSharedPtr CreateStaticCondMatrix(const MatrixKey &mkey)
Definition: TriExp.cpp:1335
void GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:821
Principle Orthogonal Functions .
Definition: BasisType.h:46
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Defines a specification for a set of points.
Definition: Points.h:58
double NekDouble
Expansion3DSharedPtr GetLeftAdjacentElementExp() const
Definition: Expansion2D.h:223
virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const
Definition: TriExp.cpp:578
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
Definition: StdExpansion.h:165
virtual const LibUtilities::BasisSharedPtr & v_GetBasis(int dir) const
Definition: TriExp.cpp:1021
virtual void v_GeneralMatrixOp_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1526
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
boost::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
Definition: NekTypeDefs.hpp:72
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition: Expansion.cpp:85
virtual void v_FwdTrans_BndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:263
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: TriExp.cpp:625
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:1730
virtual StdRegions::Orientation v_GetCartesianEorient(int edge)
Definition: TriExp.cpp:1015
SpatialDomains::GeometrySharedPtr GetGeom() const
Definition: Expansion.cpp:150
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
void GetEdgeToElementMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, int P=-1)
Definition: StdExpansion.h:846
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
virtual DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey)
Definition: TriExp.cpp:1450
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:213
Geometry is straight-sided with constant geometric factors.
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
void HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
unsigned int GetNumPoints() const
Definition: Points.h:106
virtual void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1493
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1792
virtual void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1468
boost::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
virtual void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords)
Definition: TriExp.cpp:586
virtual StdRegions::Orientation v_GetEorient(int edge)
Definition: TriExp.cpp:1009
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
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Definition: TriExp.cpp:394
virtual void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1476
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)
GeomType
Indicates the type of element geometry.
boost::shared_ptr< Basis > BasisSharedPtr
virtual void v_LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp)
Definition: TriExp.cpp:1548
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
void v_DropLocStaticCondMatrix(const MatrixKey &mkey)
Definition: TriExp.cpp:1461
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:216
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
virtual void v_GetEdgePhysMap(const int edge, Array< OneD, int > &outarray)
Definition: TriExp.cpp:764
Array< OneD, LibUtilities::BasisSharedPtr > m_base
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Geometry is curved or has non-constant factors.
void GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:816
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1057
Describes the specification for a Basis.
Definition: Basis.h:50
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:50
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
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
virtual void v_GetEdgeInterpVals(const int edge, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:746
virtual DNekScalMatSharedPtr CreateMatrix(const MatrixKey &mkey)
Definition: TriExp.cpp:1068