Nektar++
StdTriExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File StdTriExp.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: Triangle routines built upon StdExpansion2D
33 //
34 ///////////////////////////////////////////////////////////////////////////////
36 #include <StdRegions/StdTriExp.h>
38 #include <StdRegions/StdSegExp.h> // for StdSegExp, etc
39 
40 namespace Nektar
41 {
42  namespace StdRegions
43  {
44 
46  {
47  }
48 
49 
51  const LibUtilities::BasisKey &Ba,
52  const LibUtilities::BasisKey &Bb) :
53  StdExpansion (LibUtilities::StdTriData::getNumberOfCoefficients(
54  Ba.GetNumModes(),
55  Bb.GetNumModes()),
56  2,Ba,Bb),
57  StdExpansion2D(LibUtilities::StdTriData::getNumberOfCoefficients(
58  Ba.GetNumModes(),
59  Bb.GetNumModes()),
60  Ba,Bb)
61  {
62  ASSERTL0(Ba.GetNumModes() <= Bb.GetNumModes(),
63  "order in 'a' direction is higher than order "
64  "in 'b' direction");
65  }
66 
68  StdExpansion(T),
70  {
71  }
72 
74  {
75  }
76 
77  //-------------------------------
78  // Integration Methods
79  //-------------------------------
81  const Array<OneD, const NekDouble>& inarray)
82  {
83  int i;
84  int nquad1 = m_base[1]->GetNumPoints();
85  Array<OneD, NekDouble> w1_tmp(nquad1);
86 
87  Array<OneD, const NekDouble> w0 = m_base[0]->GetW();
88  Array<OneD, const NekDouble> w1 = m_base[1]->GetW();
89  Array<OneD, const NekDouble> z1 = m_base[1]->GetZ();
90 
91  switch(m_base[1]->GetPointsType())
92  {
93  case LibUtilities::eGaussRadauMAlpha1Beta0: // (0,1) Jacobi Inner product
94  {
95  Vmath::Smul(nquad1, 0.5, w1, 1, w1_tmp,1);
96  break;
97  }
98  default:
99  {
100  // include jacobian factor on whatever coordinates are defined.
101  for(i = 0; i < nquad1; ++i)
102  {
103  w1_tmp[i] = 0.5*(1-z1[i])*w1[i];
104  }
105  break;
106  }
107  }
108 
109  return StdExpansion2D::Integral(inarray,w0,w1_tmp);
110  }
111 
112  //-----------------------------
113  // Differentiation Methods
114  //-----------------------------
115 
116  /**
117  * \brief Calculate the derivative of the physical points.
118  *
119  * \f$ \frac{\partial u}{\partial x_1} = \left .
120  * \frac{2.0}{1-\eta_2} \frac{\partial u}{\partial d\eta_1}
121  * \right |_{\eta_2}\f$
122  *
123  * \f$ \frac{\partial u}{\partial x_2} = \left .
124  * \frac{1+\eta_1}{1-\eta_2} \frac{\partial u}{\partial d\eta_1}
125  * \right |_{\eta_2} + \left . \frac{\partial u}{\partial d\eta_2}
126  * \right |_{\eta_1} \f$
127  */
129  const Array<OneD, const NekDouble>& inarray,
130  Array<OneD, NekDouble>& out_d0,
131  Array<OneD, NekDouble>& out_d1,
132  Array<OneD, NekDouble>& out_d2)
133  {
134  int i;
135  int nquad0 = m_base[0]->GetNumPoints();
136  int nquad1 = m_base[1]->GetNumPoints();
137  Array<OneD, NekDouble> wsp(nquad0*nquad1);
138 
139  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
140  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
141 
142  // set up geometric factor: 2/(1-z1)
143  for (i = 0; i < nquad1; ++i)
144  {
145  wsp[i] = 2.0/(1-z1[i]);
146  }
147 
148  if (out_d0.num_elements() > 0)
149  {
150  PhysTensorDeriv(inarray, out_d0, out_d1);
151 
152  for (i = 0; i < nquad1; ++i)
153  {
154  Blas::Dscal(nquad0,wsp[i],&out_d0[0]+i*nquad0,1);
155  }
156 
157  // if no d1 required do not need to calculate both deriv
158  if (out_d1.num_elements() > 0)
159  {
160  // set up geometric factor: (1_z0)/(1-z1)
161  for (i = 0; i < nquad0; ++i)
162  {
163  wsp[i] = 0.5*(1+z0[i]);
164  }
165 
166  for (i = 0; i < nquad1; ++i)
167  {
168  Vmath::Vvtvp(nquad0,&wsp[0],1,&out_d0[0]+i*nquad0,
169  1,&out_d1[0]+i*nquad0,
170  1,&out_d1[0]+i*nquad0,1);
171  }
172  }
173  }
174  else if (out_d1.num_elements() > 0)
175  {
176  Array<OneD, NekDouble> diff0(nquad0*nquad1);
177  PhysTensorDeriv(inarray, diff0, out_d1);
178 
179  for (i = 0; i < nquad1; ++i)
180  {
181  Blas::Dscal(nquad0,wsp[i],&diff0[0]+i*nquad0,1);
182  }
183 
184  for (i = 0; i < nquad0; ++i)
185  {
186  wsp[i] = 0.5*(1+z0[i]);
187  }
188 
189  for (i = 0; i < nquad1; ++i)
190  {
191  Vmath::Vvtvp(nquad0,&wsp[0],1,&diff0[0]+i*nquad0,
192  1,&out_d1[0]+i*nquad0,
193  1,&out_d1[0]+i*nquad0,1);
194  }
195  }
196  }
197 
199  const int dir,
200  const Array<OneD, const NekDouble>& inarray,
201  Array<OneD, NekDouble>& outarray)
202  {
203  switch(dir)
204  {
205  case 0:
206  {
207  v_PhysDeriv(inarray, outarray, NullNekDouble1DArray);
208  break;
209  }
210  case 1:
211  {
212  v_PhysDeriv(inarray, NullNekDouble1DArray, outarray);
213  break;
214  }
215  default:
216  {
217  ASSERTL1(false,"input dir is out of range");
218  break;
219  }
220  }
221  }
222 
224  const Array<OneD, const NekDouble>& inarray,
225  Array<OneD, NekDouble>& out_d0,
226  Array<OneD, NekDouble>& out_d1,
227  Array<OneD, NekDouble>& out_d2)
228  {
229  StdTriExp::v_PhysDeriv(inarray, out_d0, out_d1);
230  }
231 
233  const int dir,
234  const Array<OneD, const NekDouble>& inarray,
235  Array<OneD, NekDouble>& outarray)
236  {
237  StdTriExp::v_PhysDeriv(dir,inarray,outarray);
238  }
239 
240 
241  //---------------------------------------
242  // Transforms
243  //---------------------------------------
244 
245  /**
246  * \brief Backward tranform for triangular elements
247  *
248  * @note 'q' (base[1]) runs fastest in this element.
249  */
251  const Array<OneD, const NekDouble>& inarray,
252  Array<OneD, NekDouble>& outarray)
253  {
254  v_BwdTrans_SumFac(inarray,outarray);
255  }
256 
257 
259  const Array<OneD, const NekDouble>& inarray,
260  Array<OneD, NekDouble>& outarray)
261  {
263  m_base[1]->GetNumModes());
264 
265  BwdTrans_SumFacKernel(m_base[0]->GetBdata(),
266  m_base[1]->GetBdata(),
267  inarray,outarray,wsp);
268  }
269 
271  const Array<OneD, const NekDouble>& base0,
272  const Array<OneD, const NekDouble>& base1,
273  const Array<OneD, const NekDouble>& inarray,
274  Array<OneD, NekDouble>& outarray,
276  bool doCheckCollDir0,
277  bool doCheckCollDir1)
278  {
279  int i;
280  int mode;
281  int nquad0 = m_base[0]->GetNumPoints();
282  int nquad1 = m_base[1]->GetNumPoints();
283  int nmodes0 = m_base[0]->GetNumModes();
284  int nmodes1 = m_base[1]->GetNumModes();
285 
286  ASSERTL1(wsp.num_elements() >= nquad0*nmodes1,
287  "Workspace size is not sufficient");
290  "Basis[1] is not of general tensor type");
291 
292  for (i = mode = 0; i < nmodes0; ++i)
293  {
294  Blas::Dgemv('N', nquad1,nmodes1-i,1.0,base1.get()+mode*nquad1,
295  nquad1,&inarray[0]+mode,1,0.0,&wsp[0]+i*nquad1,1);
296  mode += nmodes1-i;
297  }
298 
299  // fix for modified basis by splitting top vertex mode
301  {
302  Blas::Daxpy(nquad1,inarray[1],base1.get()+nquad1,1,
303  &wsp[0]+nquad1,1);
304  }
305 
306  Blas::Dgemm('N','T', nquad0,nquad1,nmodes0,1.0, base0.get(),nquad0,
307  &wsp[0], nquad1,0.0, &outarray[0], nquad0);
308  }
309 
311  const Array<OneD, const NekDouble>& inarray,
312  Array<OneD, NekDouble>& outarray)
313  {
314  v_IProductWRTBase(inarray,outarray);
315 
316  // get Mass matrix inverse
317  StdMatrixKey masskey(eInvMass,DetShapeType(),*this);
318  DNekMatSharedPtr matsys = GetStdMatrix(masskey);
319 
320  // copy inarray in case inarray == outarray
321  NekVector<NekDouble> in(m_ncoeffs,outarray,eCopy);
323 
324  out = (*matsys)*in;
325  }
326 
327 
329  const Array<OneD, const NekDouble>& inarray,
330  Array<OneD, NekDouble>& outarray)
331  {
332  int i,j;
333  int npoints[2] = {m_base[0]->GetNumPoints(),
334  m_base[1]->GetNumPoints()};
335  int nmodes[2] = {m_base[0]->GetNumModes(),
336  m_base[1]->GetNumModes()};
337 
338  fill(outarray.get(), outarray.get()+m_ncoeffs, 0.0 );
339 
340  Array<OneD, NekDouble> physEdge[3];
341  Array<OneD, NekDouble> coeffEdge[3];
342  for(i = 0; i < 3; i++)
343  {
344  physEdge[i] = Array<OneD, NekDouble>(npoints[i!=0]);
345  coeffEdge[i] = Array<OneD, NekDouble>(nmodes[i!=0]);
346  }
347 
348  for(i = 0; i < npoints[0]; i++)
349  {
350  physEdge[0][i] = inarray[i];
351  }
352 
353  for(i = 0; i < npoints[1]; i++)
354  {
355  physEdge[1][i] = inarray[npoints[0]-1+i*npoints[0]];
356  physEdge[2][i] = inarray[(npoints[1]-1)*npoints[0]-i*npoints[0]];
357  }
358 
359  StdSegExpSharedPtr segexp[2] = {
361  m_base[0]->GetBasisKey()),
363  m_base[1]->GetBasisKey())
364  };
365 
366  Array<OneD, unsigned int> mapArray;
367  Array<OneD, int> signArray;
368  NekDouble sign;
369 
370  for (i = 0; i < 3; i++)
371  {
372  //segexp[i!=0]->v_FwdTrans_BndConstrained(physEdge[i],coeffEdge[i]);
373  segexp[i!=0]->FwdTrans_BndConstrained(physEdge[i],coeffEdge[i]);
374 
375  v_GetEdgeToElementMap(i,eForwards,mapArray,signArray);
376  for (j = 0; j < nmodes[i != 0]; j++)
377  {
378  sign = (NekDouble) signArray[j];
379  outarray[ mapArray[j] ] = sign * coeffEdge[i][j];
380  }
381  }
382 
385 
386  StdMatrixKey masskey(eMass,DetShapeType(),*this);
387  MassMatrixOp(outarray,tmp0,masskey);
388  v_IProductWRTBase(inarray,tmp1);
389 
390  Vmath::Vsub(m_ncoeffs, tmp1, 1, tmp0, 1, tmp1, 1);
391 
392  // get Mass matrix inverse (only of interior DOF)
393  // use block (1,1) of the static condensed system
394  // note: this block alreay contains the inverse matrix
395  DNekMatSharedPtr matsys =
396  (m_stdStaticCondMatrixManager[masskey])->GetBlock(1,1);
397 
398  int nBoundaryDofs = v_NumBndryCoeffs();
399  int nInteriorDofs = m_ncoeffs - nBoundaryDofs;
400 
401  Array<OneD, NekDouble> rhs (nInteriorDofs);
402  Array<OneD, NekDouble> result(nInteriorDofs);
403 
404  v_GetInteriorMap(mapArray);
405 
406  for (i = 0; i < nInteriorDofs; i++)
407  {
408  rhs[i] = tmp1[ mapArray[i] ];
409  }
410 
411  Blas::Dgemv('N',nInteriorDofs,nInteriorDofs,
412  1.0,&(matsys->GetPtr())[0],nInteriorDofs,
413  rhs.get(),1,
414  0.0,result.get(),1);
415 
416  for (i = 0; i < nInteriorDofs; i++)
417  {
418  outarray[ mapArray[i] ] = result[i];
419  }
420  }
421 
422  //---------------------------------------
423  // Inner product functions
424  //---------------------------------------
425 
426  /**
427  * \brief Calculate the inner product of inarray with respect to the
428  * basis B=base0[p]*base1[pq] and put into outarray.
429  *
430  * \f$
431  * \begin{array}{rcl}
432  * I_{pq} = (\phi^A_q \phi^B_{pq}, u) &=&
433  * \sum_{i=0}^{nq_0}\sum_{j=0}^{nq_1}
434  * \phi^A_p(\eta_{0,i})\phi^B_{pq}(\eta_{1,j}) w^0_i w^1_j
435  * u(\xi_{0,i} \xi_{1,j}) \\
436  * & = & \sum_{i=0}^{nq_0} \phi^A_p(\eta_{0,i})
437  * \sum_{j=0}^{nq_1} \phi^B_{pq}(\eta_{1,j}) \tilde{u}_{i,j}
438  * \end{array}
439  * \f$
440  *
441  * where
442  *
443  * \f$ \tilde{u}_{i,j} = w^0_i w^1_j u(\xi_{0,i},\xi_{1,j}) \f$
444  *
445  * which can be implemented as
446  *
447  * \f$ f_{pj} = \sum_{i=0}^{nq_0} \phi^A_p(\eta_{0,i})
448  * \tilde{u}_{i,j}
449  * \rightarrow {\bf B_1 U} \f$
450  * \f$ I_{pq} = \sum_{j=0}^{nq_1} \phi^B_{pq}(\eta_{1,j}) f_{pj}
451  * \rightarrow {\bf B_2[p*skip] f[skip]} \f$
452  *
453  * \b Recall: \f$ \eta_{1} = \frac{2(1+\xi_1)}{(1-\xi_2)}-1, \,
454  * \eta_2 = \xi_2\f$
455  *
456  * \b Note: For the orthgonality of this expansion to be realised the
457  * 'q' ordering must run fastest in contrast to the Quad and Hex
458  * ordering where 'p' index runs fastest to be consistent with the
459  * quadrature ordering.
460  *
461  * In the triangular space the i (i.e. \f$\eta_1\f$ direction)
462  * ordering still runs fastest by convention.
463  */
465  const Array<OneD, const NekDouble>& inarray,
466  Array<OneD, NekDouble>& outarray)
467  {
468  StdTriExp::v_IProductWRTBase_SumFac(inarray,outarray);
469  }
470 
472  const Array<OneD, const NekDouble>& inarray,
473  Array<OneD, NekDouble>& outarray)
474  {
475  int nq = GetTotPoints();
476  StdMatrixKey iprodmatkey(eIProductWRTBase,DetShapeType(),*this);
477  DNekMatSharedPtr iprodmat = GetStdMatrix(iprodmatkey);
478 
479  Blas::Dgemv('N',m_ncoeffs,nq,1.0,iprodmat->GetPtr().get(),
480  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
481  }
482 
484  const Array<OneD, const NekDouble>& inarray,
485  Array<OneD, NekDouble>& outarray,
486  bool multiplybyweights)
487  {
488  int nquad0 = m_base[0]->GetNumPoints();
489  int nquad1 = m_base[1]->GetNumPoints();
490  int order0 = m_base[0]->GetNumModes();
491 
492  if(multiplybyweights)
493  {
494  Array<OneD,NekDouble> tmp(nquad0*nquad1+nquad1*order0);
495  Array<OneD,NekDouble> wsp(tmp+nquad0*nquad1);
496 
497  // multiply by integration constants
498  MultiplyByQuadratureMetric(inarray,tmp);
499  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
500  m_base[1]->GetBdata(),
501  tmp,outarray,wsp);
502  }
503  else
504  {
505  Array<OneD,NekDouble> wsp(nquad1*order0);
506  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
507  m_base[1]->GetBdata(),
508  inarray,outarray,wsp);
509  }
510  }
511 
513  const Array<OneD, const NekDouble>& base0,
514  const Array<OneD, const NekDouble>& base1,
515  const Array<OneD, const NekDouble>& inarray,
516  Array<OneD, NekDouble>& outarray,
518  bool doCheckCollDir0,
519  bool doCheckCollDir1)
520  {
521  int i;
522  int mode;
523  int nquad0 = m_base[0]->GetNumPoints();
524  int nquad1 = m_base[1]->GetNumPoints();
525  int nmodes0 = m_base[0]->GetNumModes();
526  int nmodes1 = m_base[1]->GetNumModes();
527 
528  ASSERTL1(wsp.num_elements() >= nquad1*nmodes0,
529  "Workspace size is not sufficient");
530 
531  Blas::Dgemm('T','N',nquad1,nmodes0,nquad0,1.0,inarray.get(),nquad0,
532  base0.get(),nquad0,0.0,wsp.get(),nquad1);
533 
534  // Inner product with respect to 'b' direction
535  for (mode=i=0; i < nmodes0; ++i)
536  {
537  Blas::Dgemv('T',nquad1,nmodes1-i,1.0, base1.get()+mode*nquad1,
538  nquad1,wsp.get()+i*nquad1,1, 0.0,
539  outarray.get() + mode,1);
540  mode += nmodes1 - i;
541  }
542 
543  // fix for modified basis by splitting top vertex mode
545  {
546  outarray[1] += Blas::Ddot(nquad1,base1.get()+nquad1,1,
547  wsp.get()+nquad1,1);
548  }
549  }
550 
552  const int dir,
553  const Array<OneD, const NekDouble>& inarray,
554  Array<OneD, NekDouble>& outarray)
555  {
556  StdTriExp::v_IProductWRTDerivBase_SumFac(dir,inarray,outarray);
557  }
558 
560  const int dir,
561  const Array<OneD, const NekDouble>& inarray,
562  Array<OneD, NekDouble>& outarray)
563  {
564  int nq = GetTotPoints();
566 
567  switch(dir)
568  {
569  case 0:
570  {
571  mtype = eIProductWRTDerivBase0;
572  break;
573  }
574  case 1:
575  {
576  mtype = eIProductWRTDerivBase1;
577  break;
578  }
579  default:
580  {
581  ASSERTL1(false,"input dir is out of range");
582  break;
583  }
584  }
585 
586  StdMatrixKey iprodmatkey(mtype,DetShapeType(),*this);
587  DNekMatSharedPtr iprodmat = GetStdMatrix(iprodmatkey);
588 
589  Blas::Dgemv('N',m_ncoeffs,nq,1.0,iprodmat->GetPtr().get(),
590  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
591  }
592 
594  const int dir,
595  const Array<OneD, const NekDouble>& inarray,
596  Array<OneD, NekDouble>& outarray)
597  {
598  int i;
599  int nquad0 = m_base[0]->GetNumPoints();
600  int nquad1 = m_base[1]->GetNumPoints();
601  int nqtot = nquad0*nquad1;
602  int nmodes0 = m_base[0]->GetNumModes();
603  int wspsize = max(max(nqtot,m_ncoeffs),nquad1*nmodes0);
604 
605  Array<OneD, NekDouble> gfac0(2*wspsize);
606  Array<OneD, NekDouble> tmp0 (gfac0+wspsize);
607 
608  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
609 
610  // set up geometric factor: 2/(1-z1)
611  for(i = 0; i < nquad1; ++i)
612  {
613  gfac0[i] = 2.0/(1-z1[i]);
614  }
615 
616  for(i = 0; i < nquad1; ++i)
617  {
618  Vmath::Smul(nquad0,gfac0[i],&inarray[0]+i*nquad0,1,
619  &tmp0[0]+i*nquad0,1);
620  }
621 
622  MultiplyByQuadratureMetric(tmp0,tmp0);
623 
624  switch(dir)
625  {
626  case 0:
627  {
628  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
629  m_base[1]->GetBdata(),
630  tmp0,outarray,gfac0);
631  break;
632  }
633  case 1:
634  {
636  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
637 
638  for (i = 0; i < nquad0; ++i)
639  {
640  gfac0[i] = 0.5*(1+z0[i]);
641  }
642 
643  for (i = 0; i < nquad1; ++i)
644  {
645  Vmath::Vmul(nquad0,&gfac0[0],1,&tmp0[0]+i*nquad0,1,
646  &tmp0[0]+i*nquad0,1);
647  }
648 
649  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
650  m_base[1]->GetBdata(),
651  tmp0,tmp3,gfac0);
652 
653  MultiplyByQuadratureMetric(inarray,tmp0);
654  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
655  m_base[1]->GetDbdata(),
656  tmp0,outarray,gfac0);
657  Vmath::Vadd(m_ncoeffs,&tmp3[0],1,&outarray[0],1,
658  &outarray[0],1);
659  break;
660  }
661  default:
662  {
663  ASSERTL1(false, "input dir is out of range");
664  break;
665  }
666  }
667  }
668 
669  //---------------------------------------
670  // Evaluation functions
671  //---------------------------------------
672 
675  {
676 
677  // set up local coordinate system
678  if (fabs(xi[1]-1.0) < NekConstants::kNekZeroTol)
679  {
680  eta[0] = -1.0;
681  eta[1] = 1.0;
682  }
683  else
684  {
685  eta[0] = 2*(1+xi[0])/(1-xi[1])-1.0;
686  eta[1] = xi[1];
687  }
688  }
689 
691  const int mode, Array<OneD, NekDouble> &outarray)
692  {
693  int i,m;
694  int nquad0 = m_base[0]->GetNumPoints();
695  int nquad1 = m_base[1]->GetNumPoints();
696  int order0 = m_base[0]->GetNumModes();
697  int order1 = m_base[1]->GetNumModes();
698  int mode0 = 0;
699  Array<OneD, const NekDouble> base0 = m_base[0]->GetBdata();
700  Array<OneD, const NekDouble> base1 = m_base[1]->GetBdata();
701 
702  ASSERTL2(mode <= m_ncoeffs,
703  "calling argument mode is larger than "
704  "total expansion order");
705 
706  m = order1;
707  for (i = 0; i < order0; ++i, m+=order1-i)
708  {
709  if (m > mode)
710  {
711  mode0 = i;
712  break;
713  }
714  }
715 
716  // deal with top vertex mode in modified basis
717  if (mode == 1 &&
719  {
720  Vmath::Fill(nquad0*nquad1 , 1.0, outarray, 1);
721  }
722  else
723  {
724  for (i = 0; i < nquad1; ++i)
725  {
726  Vmath::Vcopy(nquad0,(NekDouble *)(base0.get()+mode0*nquad0),
727  1,&outarray[0]+i*nquad0,1);
728  }
729  }
730 
731  for(i = 0; i < nquad0; ++i)
732  {
733  Vmath::Vmul(nquad1,(NekDouble *)(base1.get() + mode*nquad1),
734  1,&outarray[0]+i,nquad0,&outarray[0]+i,nquad0);
735  }
736  }
737 
738 
740  {
741  return 3;
742  }
743 
745  {
746  return 3;
747  }
748 
750  {
752  }
753 
755  {
757  "BasisType is not a boundary interior form");
759  "BasisType is not a boundary interior form");
760 
761  return 3 + (GetBasisNumModes(0)-2) + 2*(GetBasisNumModes(1)-2);
762  }
763 
765  {
767  "BasisType is not a boundary interior form");
769  "BasisType is not a boundary interior form");
770 
771  return GetBasisNumModes(0) + 2*GetBasisNumModes(1);
772  }
773 
774  int StdTriExp::v_GetEdgeNcoeffs(const int i) const
775  {
776  ASSERTL2(i >= 0 && i <= 2, "edge id is out of range");
777 
778  if (i == 0)
779  {
780  return GetBasisNumModes(0);
781  }
782  else
783  {
784  return GetBasisNumModes(1);
785  }
786  }
787 
788  int StdTriExp::v_GetEdgeNumPoints(const int i) const
789  {
790  ASSERTL2((i >= 0)&&(i <= 2),"edge id is out of range");
791 
792  if (i == 0)
793  {
794  return GetNumPoints(0);
795  }
796  else
797  {
798  return GetNumPoints(1);
799  }
800  }
801 
803  const std::vector<unsigned int> &nummodes,
804  int &modes_offset)
805  {
807  nummodes[modes_offset],
808  nummodes[modes_offset+1]);
809  modes_offset += 2;
810 
811  return nmodes;
812  }
813 
815  {
816  ASSERTL2(i >= 0 && i <= 2, "edge id is out of range");
817 
818  if (i == 0)
819  {
820  return GetBasisType(0);
821  }
822  else
823  {
824  return GetBasisType(1);
825  }
826  }
827 
828 
830  Array<OneD, NekDouble> &coords_1,
831  Array<OneD, NekDouble> &coords_2)
832  {
833  Array<OneD, const NekDouble> z0 = m_base[0]->GetZ();
834  Array<OneD, const NekDouble> z1 = m_base[1]->GetZ();
835  int nq0 = GetNumPoints(0);
836  int nq1 = GetNumPoints(1);
837  int i,j;
838 
839  for(i = 0; i < nq1; ++i)
840  {
841  for(j = 0; j < nq0; ++j)
842  {
843  coords_0[i*nq0+j] = (1+z0[j])*(1-z1[i])/2.0 - 1.0;
844  }
845  Vmath::Fill(nq0,z1[i],&coords_1[0] + i*nq0,1);
846  }
847  }
848 
850  {
851  return m_base[0]->GetBasisType() == LibUtilities::eModified_A &&
852  m_base[1]->GetBasisType() == LibUtilities::eModified_B;
853  }
854 
856  {
857  ASSERTL2(edge >= 0 && edge <= 2, "edge id is out of range");
858 
859  return edge == 0 ? 0 : 1;
860  }
861 
863  const int i) const
864  {
865  ASSERTL2(i >= 0 && i <= 2, "edge id is out of range");
866 
867  if (i == 0)
868  {
869  return GetBasis(0)->GetBasisKey();
870  }
871  else
872  {
873  switch(m_base[1]->GetBasisType())
874  {
876  {
877  switch(m_base[1]->GetPointsType())
878  {
880  {
882  m_base[1]->GetBasisKey().GetPointsKey().
883  GetNumPoints()+1,
885  return LibUtilities::BasisKey(
887  m_base[1]->GetNumModes(),pkey);
888  break;
889  }
890 
891  default:
892  ASSERTL0(false,"unexpected points distribution");
893  break;
894  }
895  }
896  default:
897  ASSERTL0(false,"Information not available to set edge key");
898  break;
899  }
900  }
902  }
903 
904 
905 
906  //--------------------------
907  // Mappings
908  //--------------------------
909 
911  const int eid,
912  const Orientation edgeOrient,
913  Array<OneD, unsigned int>& maparray,
914  Array<OneD, int>& signarray)
915  {
918  "Mapping not defined for this type of basis");
919 
920  int i;
921  const int nummodes1 = m_base[1]->GetNumModes();
922  const int nEdgeCoeffs = GetEdgeNcoeffs(eid);
923 
924  if(maparray.num_elements() != nEdgeCoeffs)
925  {
926  maparray = Array<OneD, unsigned int>(nEdgeCoeffs);
927  }
928 
929  if(signarray.num_elements() != nEdgeCoeffs)
930  {
931  signarray = Array<OneD, int>(nEdgeCoeffs,1);
932  }
933  else
934  {
935  fill(signarray.get() , signarray.get()+nEdgeCoeffs, 1);
936  }
937 
938  switch(eid)
939  {
940  case 0:
941  {
942  int cnt = 0;
943  for(i = 0; i < nEdgeCoeffs; cnt+=nummodes1-i, ++i)
944  {
945  maparray[i] = cnt;
946  }
947 
948  if(edgeOrient==eBackwards)
949  {
950  swap( maparray[0] , maparray[1] );
951 
952  for(i = 3; i < nEdgeCoeffs; i+=2)
953  {
954  signarray[i] = -1;
955  }
956  }
957  break;
958  }
959  case 1:
960  {
961  maparray[0] = nummodes1;
962  maparray[1] = 1;
963  for(i = 2; i < nEdgeCoeffs; i++)
964  {
965  maparray[i] = nummodes1-1+i;
966  }
967 
968  if(edgeOrient==eBackwards)
969  {
970  swap( maparray[0] , maparray[1] );
971 
972  for(i = 3; i < nEdgeCoeffs; i+=2)
973  {
974  signarray[i] = -1;
975  }
976  }
977  break;
978  }
979  case 2:
980  {
981  for(i = 0; i < nEdgeCoeffs; i++)
982  {
983  maparray[i] = i;
984  }
985 
986  if(edgeOrient==eForwards)
987  {
988  swap( maparray[0] , maparray[1] );
989 
990  for(i = 3; i < nEdgeCoeffs; i+=2)
991  {
992  signarray[i] = -1;
993  }
994  }
995  break;
996  }
997  default:
998  ASSERTL0(false,"eid must be between 0 and 2");
999  break;
1000  }
1001  }
1002 
1003  int StdTriExp::v_GetVertexMap(const int localVertexId,bool useCoeffPacking)
1004  {
1005  ASSERTL0(
1006  GetEdgeBasisType(localVertexId) == LibUtilities::eModified_A ||
1007  GetEdgeBasisType(localVertexId) == LibUtilities::eModified_B,
1008  "Mapping not defined for this type of basis");
1009 
1010  int localDOF = 0;
1011  if(useCoeffPacking == true)
1012  {
1013  switch(localVertexId)
1014  {
1015  case 0:
1016  {
1017  localDOF = 0;
1018  break;
1019  }
1020  case 1:
1021  {
1022  localDOF = 1;
1023  break;
1024  }
1025  case 2:
1026  {
1027  localDOF = m_base[1]->GetNumModes();
1028  break;
1029  }
1030  default:
1031  {
1032  ASSERTL0(false,"eid must be between 0 and 2");
1033  break;
1034  }
1035  }
1036  }
1037  else // follow book format for vertex indexing.
1038  {
1039  switch(localVertexId)
1040  {
1041  case 0:
1042  {
1043  localDOF = 0;
1044  break;
1045  }
1046  case 1:
1047  {
1048  localDOF = m_base[1]->GetNumModes();
1049  break;
1050  }
1051  case 2:
1052  {
1053  localDOF = 1;
1054  break;
1055  }
1056  default:
1057  {
1058  ASSERTL0(false,"eid must be between 0 and 2");
1059  break;
1060  }
1061  }
1062  }
1063 
1064  return localDOF;
1065  }
1066 
1068  const int eid,
1069  const Orientation edgeOrient,
1070  Array<OneD, unsigned int>& maparray,
1071  Array<OneD, int>& signarray)
1072  {
1075  "Mapping not defined for this type of basis");
1076  int i;
1077  const int nummodes1 = m_base[1]->GetNumModes();
1078  const int nEdgeIntCoeffs = GetEdgeNcoeffs(eid)-2;
1079 
1080  if(maparray.num_elements() != nEdgeIntCoeffs)
1081  {
1082  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1083  }
1084 
1085  if(signarray.num_elements() != nEdgeIntCoeffs)
1086  {
1087  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
1088  }
1089  else
1090  {
1091  fill( signarray.get() , signarray.get()+nEdgeIntCoeffs, 1 );
1092  }
1093 
1094  switch(eid)
1095  {
1096  case 0:
1097  {
1098  int cnt = 2*nummodes1 - 1;
1099  for(i = 0; i < nEdgeIntCoeffs; cnt+=nummodes1-2-i, ++i)
1100  {
1101  maparray[i] = cnt;
1102  }
1103 
1104  if(edgeOrient==eBackwards)
1105  {
1106  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1107  {
1108  signarray[i] = -1;
1109  }
1110  }
1111  break;
1112  }
1113  case 1:
1114  {
1115  for(i = 0; i < nEdgeIntCoeffs; i++)
1116  {
1117  maparray[i] = nummodes1+1+i;
1118  }
1119 
1120  if(edgeOrient==eBackwards)
1121  {
1122  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1123  {
1124  signarray[i] = -1;
1125  }
1126  }
1127  break;
1128  }
1129  case 2:
1130  {
1131  for(i = 0; i < nEdgeIntCoeffs; i++)
1132  {
1133  maparray[i] = 2+i;
1134  }
1135 
1136  if(edgeOrient==eForwards)
1137  {
1138  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1139  {
1140  signarray[i] = -1;
1141  }
1142  }
1143  break;
1144  }
1145  default:
1146  {
1147  ASSERTL0(false,"eid must be between 0 and 2");
1148  break;
1149  }
1150  }
1151  }
1152 
1154  {
1157  "Expansion not of a proper type");
1158 
1159  int i,j;
1160  int cnt = 0;
1161  int nummodes0, nummodes1;
1162  int startvalue;
1163  if(outarray.num_elements()!=GetNcoeffs()-NumBndryCoeffs())
1164  {
1166  }
1167 
1168  nummodes0 = m_base[0]->GetNumModes();
1169  nummodes1 = m_base[1]->GetNumModes();
1170 
1171  startvalue = 2*nummodes1;
1172 
1173  for(i = 0; i < nummodes0-2; i++)
1174  {
1175  for(j = 0; j < nummodes1-3-i; j++)
1176  {
1177  outarray[cnt++]=startvalue+j;
1178  }
1179  startvalue+=nummodes1-2-i;
1180  }
1181  }
1182 
1184  {
1187  "Expansion not of expected type");
1188  int i;
1189  int cnt;
1190  int nummodes0, nummodes1;
1191  int value;
1192 
1193  if (outarray.num_elements()!=NumBndryCoeffs())
1194  {
1196  }
1197 
1198  nummodes0 = m_base[0]->GetNumModes();
1199  nummodes1 = m_base[1]->GetNumModes();
1200 
1201  value = 2*nummodes1-1;
1202  for(i = 0; i < value; i++)
1203  {
1204  outarray[i]=i;
1205  }
1206  cnt = value;
1207 
1208  for(i = 0; i < nummodes0-2; i++)
1209  {
1210  outarray[cnt++]=value;
1211  value += nummodes1-2-i;
1212  }
1213  }
1214 
1215 
1216  //---------------------------------------
1217  // Wrapper functions
1218  //---------------------------------------
1219 
1221  {
1222 
1223  MatrixType mtype = mkey.GetMatrixType();
1224 
1225  DNekMatSharedPtr Mat;
1226 
1227  switch(mtype)
1228  {
1230  {
1231  int nq0 = m_base[0]->GetNumPoints();
1232  int nq1 = m_base[1]->GetNumPoints();
1233  int nq = max(nq0,nq1);
1234  int neq = LibUtilities::StdTriData::
1236  Array<OneD, Array<OneD, NekDouble> > coords(neq);
1237  Array<OneD, NekDouble> coll (2);
1239  Array<OneD, NekDouble> tmp (nq0);
1240 
1241  Mat = MemoryManager<DNekMat>::AllocateSharedPtr(neq,nq0*nq1);
1242  int cnt = 0;
1243 
1244  for(int i = 0; i < nq; ++i)
1245  {
1246  for(int j = 0; j < nq-i; ++j,++cnt)
1247  {
1248  coords[cnt] = Array<OneD, NekDouble>(2);
1249  coords[cnt][0] = -1.0 + 2*j/(NekDouble)(nq-1);
1250  coords[cnt][1] = -1.0 + 2*i/(NekDouble)(nq-1);
1251  }
1252  }
1253 
1254  for(int i = 0; i < neq; ++i)
1255  {
1256  LocCoordToLocCollapsed(coords[i],coll);
1257 
1258  I[0] = m_base[0]->GetI(coll);
1259  I[1] = m_base[1]->GetI(coll+1);
1260 
1261  // interpolate first coordinate direction
1262  for (int j = 0; j < nq1; ++j)
1263  {
1264  NekDouble fac = (I[1]->GetPtr())[j];
1265  Vmath::Smul(nq0, fac, I[0]->GetPtr(), 1, tmp, 1);
1266 
1267  Vmath::Vcopy(nq0, &tmp[0], 1,
1268  Mat->GetRawPtr() + j*nq0*neq + i, neq);
1269  }
1270 
1271  }
1272  break;
1273  }
1274  default:
1275  {
1277  break;
1278  }
1279  }
1280 
1281  return Mat;
1282  }
1283 
1285  {
1286  return v_GenMatrix(mkey);
1287  }
1288 
1289 
1290  //---------------------------------------
1291  // Operator evaluation functions
1292  //---------------------------------------
1293 
1295  const Array<OneD, const NekDouble> &inarray,
1296  Array<OneD, NekDouble> &outarray,
1297  const StdMatrixKey &mkey)
1298  {
1299  StdExpansion::MassMatrixOp_MatFree(inarray,outarray,mkey);
1300  }
1301 
1303  const Array<OneD, const NekDouble> &inarray,
1304  Array<OneD, NekDouble> &outarray,
1305  const StdMatrixKey &mkey)
1306  {
1307  StdTriExp::v_LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
1308  }
1309 
1311  const int k1,
1312  const int k2,
1313  const Array<OneD, const NekDouble> &inarray,
1314  Array<OneD, NekDouble> &outarray,
1315  const StdMatrixKey &mkey)
1316  {
1318  k1,k2,inarray,outarray,mkey);
1319  }
1320 
1322  const int i,
1323  const Array<OneD, const NekDouble> &inarray,
1324  Array<OneD, NekDouble> &outarray,
1325  const StdMatrixKey &mkey)
1326  {
1327  StdExpansion::WeakDerivMatrixOp_MatFree(i,inarray,outarray,mkey);
1328  }
1329 
1331  const Array<OneD, const NekDouble> &inarray,
1332  Array<OneD, NekDouble> &outarray,
1333  const StdMatrixKey &mkey)
1334  {
1335  StdTriExp::v_HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
1336  }
1337 
1338 
1340  const StdMatrixKey &mkey)
1341  {
1342  int qa = m_base[0]->GetNumPoints();
1343  int qb = m_base[1]->GetNumPoints();
1344  int nmodes_a = m_base[0]->GetNumModes();
1345  int nmodes_b = m_base[1]->GetNumModes();
1346 
1347  // Declare orthogonal basis.
1350 
1353  StdTriExp OrthoExp(Ba,Bb);
1354 
1355  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
1356  int j, k , cnt = 0;
1357 
1358  int cutoff = (int) (mkey.GetConstFactor(eFactorSVVCutoffRatio)*min(nmodes_a,nmodes_b));
1359  NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDiffCoeff);
1360 
1361  NekDouble epsilon = 1.0;
1362  int nmodes = min(nmodes_a,nmodes_b);
1363 
1364  //To avoid the fac[j] from blowing up
1365  //NekDouble epsilon = 0.001;
1366 
1367  // project onto physical space.
1368  OrthoExp.FwdTrans(array,orthocoeffs);
1369 
1370  //cout << "nmodes_a = " << nmodes_a << " and nmodes_b = " << nmodes_b << "and and orthocoeffs is of size " << sizeof(orthocoeffs) << endl;
1371  // apply SVV filter (JEL)
1372  for(j = 0; j < nmodes_a; ++j)
1373  {
1374  for(k = 0; k < nmodes_b-j; ++k)
1375  {
1376  if(j + k >= cutoff)
1377  {
1378  orthocoeffs[cnt] *= (SvvDiffCoeff*exp(-(j+k-nmodes)*(j+k-nmodes)/((NekDouble)((j+k-cutoff+epsilon)*(j+k-cutoff+epsilon)))));
1379  }
1380  else
1381  {
1382  orthocoeffs[cnt] *= 0.0;
1383  }
1384  cnt++;
1385  }
1386  }
1387 
1388  // backward transform to physical space
1389  OrthoExp.BwdTrans(orthocoeffs,array);
1390  }
1391 
1393  int numMin,
1394  const Array<OneD, const NekDouble> &inarray,
1395  Array<OneD, NekDouble> &outarray)
1396  {
1397  int n_coeffs = inarray.num_elements();
1398  int nquad0 = m_base[0]->GetNumPoints();
1399  int nquad1 = m_base[1]->GetNumPoints();
1400  Array<OneD, NekDouble> coeff(n_coeffs);
1401  Array<OneD, NekDouble> coeff_tmp(n_coeffs,0.0);
1404  int nqtot = nquad0*nquad1;
1405  Array<OneD, NekDouble> phys_tmp(nqtot,0.0);
1406 
1407  int nmodes0 = m_base[0]->GetNumModes();
1408  int nmodes1 = m_base[1]->GetNumModes();
1409  int numMin2 = nmodes0;
1410  int i;
1411 
1412  const LibUtilities::PointsKey Pkey0(
1414  const LibUtilities::PointsKey Pkey1(
1416 
1417  LibUtilities::BasisKey b0(m_base[0]->GetBasisType(),nmodes0,Pkey0);
1418  LibUtilities::BasisKey b1(m_base[1]->GetBasisType(),nmodes1,Pkey1);
1419 
1420  LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A,nmodes0,Pkey0);
1421  LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_B,nmodes1,Pkey1);
1422 
1423  StdRegions::StdTriExpSharedPtr m_OrthoTriExp;
1425 
1427  ::AllocateSharedPtr(b0, b1);
1428  m_OrthoTriExp = MemoryManager<StdRegions::StdTriExp>
1429  ::AllocateSharedPtr(bortho0, bortho1);
1430 
1431  m_TriExp ->BwdTrans(inarray,phys_tmp);
1432  m_OrthoTriExp->FwdTrans(phys_tmp, coeff);
1433 
1434  for (i = 0; i < n_coeffs; i++)
1435  {
1436  if (i == numMin)
1437  {
1438  coeff[i] = 0.0;
1439  numMin += numMin2 - 1;
1440  numMin2 -= 1.0;
1441  }
1442  }
1443 
1444  m_OrthoTriExp->BwdTrans(coeff,phys_tmp);
1445  m_TriExp ->FwdTrans(phys_tmp, outarray);
1446 
1447  }
1448 
1450  const Array<OneD, const NekDouble> &inarray,
1451  Array<OneD, NekDouble> &outarray,
1452  const StdMatrixKey &mkey)
1453  {
1455 
1456  if(inarray.get() == outarray.get())
1457  {
1459  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
1460 
1461  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, 1.0, mat->GetPtr().get(),
1462  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
1463  }
1464  else
1465  {
1466  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, 1.0, mat->GetPtr().get(),
1467  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
1468  }
1469  }
1470 
1471  //---------------------------------------
1472  // Private helper functions
1473  //---------------------------------------
1474 
1476  const Array<OneD, const NekDouble>& inarray,
1477  Array<OneD, NekDouble> &outarray)
1478  {
1479  int i;
1480  int nquad0 = m_base[0]->GetNumPoints();
1481  int nquad1 = m_base[1]->GetNumPoints();
1482 
1483  const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
1484  const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
1485  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
1486 
1487  // multiply by integration constants
1488  for(i = 0; i < nquad1; ++i)
1489  {
1490  Vmath::Vmul(nquad0,inarray.get()+i*nquad0,1,
1491  w0.get(),1, outarray.get()+i*nquad0,1);
1492  }
1493 
1494  switch(m_base[1]->GetPointsType())
1495  {
1496  // Legendre inner product
1499  for(i = 0; i < nquad1; ++i)
1500  {
1501  Blas::Dscal(nquad0,0.5*(1-z1[i])*w1[i],
1502  outarray.get()+i*nquad0,1);
1503  }
1504  break;
1505 
1506  // (1,0) Jacobi Inner product
1508  for(i = 0; i < nquad1; ++i)
1509  {
1510  Blas::Dscal(nquad0,0.5*w1[i],outarray.get()+i*nquad0,1);
1511  }
1512  break;
1513 
1514  default:
1515  ASSERTL0(false, "Unsupported quadrature points type.");
1516  break;
1517  }
1518  }
1519 
1521  Array<OneD, int> &conn,
1522  bool standard)
1523  {
1524  int np1 = m_base[0]->GetNumPoints();
1525  int np2 = m_base[1]->GetNumPoints();
1526  int np = max(np1,np2);
1527 
1528  conn = Array<OneD, int>(3*(np-1)*(np-1));
1529 
1530  int row = 0;
1531  int rowp1 = 0;
1532  int cnt = 0;
1533  for(int i = 0; i < np-1; ++i)
1534  {
1535  rowp1 += np-i;
1536  for(int j = 0; j < np-i-2; ++j)
1537  {
1538  conn[cnt++] = row +j;
1539  conn[cnt++] = row +j+1;
1540  conn[cnt++] = rowp1 +j;
1541 
1542  conn[cnt++] = rowp1 +j+1;
1543  conn[cnt++] = rowp1 +j;
1544  conn[cnt++] = row +j+1;
1545  }
1546 
1547  conn[cnt++] = row +np-i-2;
1548  conn[cnt++] = row +np-i-1;
1549  conn[cnt++] = rowp1+np-i-2;
1550 
1551  row += np-i;
1552  }
1553  }
1554 
1555 
1556  }//end namespace
1557 }//end namespace
boost::shared_ptr< StdTriExp > StdTriExpSharedPtr
Definition: StdTriExp.h:266
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:454
virtual void v_IProductWRTDerivBase_MatOp(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:559
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:122
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
const LibUtilities::BasisSharedPtr & GetBasis(int dir) const
This function gets the shared point to basis in the dir direction.
Definition: StdExpansion.h:118
virtual void v_GeneralMatrixOp_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdTriExp.cpp:1449
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:82
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Backward tranform for triangular elements.
Definition: StdTriExp.cpp:250
void MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:931
static Array< OneD, NekDouble > NullNekDouble1DArray
void BwdTrans_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)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
int GetBasisNumModes(const int dir) const
This function returns the number of expansion modes in the dir direction.
Definition: StdExpansion.h:178
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:22
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:902
virtual void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:1475
virtual LibUtilities::ShapeType v_DetShapeType() const
Definition: StdTriExp.cpp:749
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
virtual int v_NumDGBndryCoeffs() const
Definition: StdTriExp.cpp:764
virtual void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdTriExp.cpp:1321
void LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
virtual void v_FwdTrans_BndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:328
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: StdTriExp.cpp:464
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
Definition: StdExpansion.h:229
virtual void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdTriExp.cpp:1330
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 LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
Convert local cartesian coordinate xi into local collapsed coordinates eta.
NekDouble Integral(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &w0, const Array< OneD, const NekDouble > &w1)
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:593
virtual void v_LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
LibUtilities::NekManager< StdMatrixKey, DNekBlkMat, StdMatrixKey::opLess > m_stdStaticCondMatrixManager
virtual int v_GetEdgeNumPoints(const int i) const
Definition: StdTriExp.cpp:788
int GetEdgeNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th edge.
Definition: StdExpansion.h:287
boost::shared_ptr< StdSegExp > StdSegExpSharedPtr
Definition: StdSegExp.h:47
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdTriExp.cpp:1153
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
Definition: StdTriExp.cpp:1339
virtual void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
Definition: StdTriExp.cpp:673
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:684
void PhysTensorDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray_d0, Array< OneD, NekDouble > &outarray_d1)
Calculate the 2D derivative in the local tensor/collapsed coordinate at the physical points...
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:57
virtual int v_DetCartesianDirOfEdge(const int edge)
Definition: StdTriExp.cpp:855
virtual int v_NumBndryCoeffs() const
Definition: StdTriExp.cpp:754
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
Definition: StdTriExp.cpp:1220
Principle Orthogonal Functions .
Definition: BasisType.h:47
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:63
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:141
virtual int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset)
Definition: StdTriExp.cpp:802
virtual void v_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true)
Definition: StdTriExp.cpp:1520
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:1392
virtual int v_GetNverts() const
Definition: StdTriExp.cpp:739
LibUtilities::BasisType GetEdgeBasisType(const int i) const
This function returns the type of expansion basis on the i-th edge.
Definition: StdExpansion.h:397
virtual LibUtilities::BasisType v_GetEdgeBasisType(const int i) const
Definition: StdTriExp.cpp:814
The base class for all shapes.
Definition: StdExpansion.h:69
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdTriExp.cpp:1183
static const NekDouble kNekZeroTol
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_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdTriExp.cpp:1294
Principle Modified Functions .
Definition: BasisType.h:50
virtual void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdTriExp.cpp:1302
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
Definition: StdTriExp.cpp:1284
void WeakDerivMatrixOp_MatFree(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:111
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
Principle Orthogonal Functions .
Definition: BasisType.h:46
Defines a specification for a set of points.
Definition: Points.h:58
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Definition: StdTriExp.cpp:483
virtual void v_IProductWRTBase_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:471
double NekDouble
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:258
T Ddot(int n, const Array< OneD, const T > &w, const int incw, const Array< OneD, const T > &x, const int incx, const Array< OneD, const int > &y, const int incy)
Definition: VmathArray.hpp:434
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
Definition: StdExpansion.h:165
virtual int v_GetEdgeNcoeffs(const int i) const
Definition: StdTriExp.cpp:774
virtual void v_BwdTrans_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, bool doCheckCollDir1)
Definition: StdTriExp.cpp:270
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
virtual void v_StdPhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
Definition: StdTriExp.cpp:223
virtual NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
Integrates the specified function over the domain.
Definition: StdTriExp.cpp:80
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: StdTriExp.cpp:310
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:551
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:213
void MassMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
virtual bool v_IsBoundaryInteriorExpansion()
Definition: StdTriExp.cpp:849
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Backward transformation from coefficient space to physical space...
Definition: StdExpansion.h:509
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
Definition: StdTriExp.cpp:1003
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:690
virtual int v_GetNedges() const
Definition: StdTriExp.cpp:744
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:131
LibUtilities::NekManager< StdMatrixKey, DNekMat, StdMatrixKey::opLess > m_stdMatrixManager
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: StdTriExp.cpp:128
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)
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:84
virtual void v_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, bool doCheckCollDir1)
Definition: StdTriExp.cpp:512
virtual void v_GetEdgeInteriorMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray)
Definition: StdTriExp.cpp:1067
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
Array< OneD, LibUtilities::BasisSharedPtr > m_base
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_x, Array< OneD, NekDouble > &coords_y, Array< OneD, NekDouble > &coords_z)
Definition: StdTriExp.cpp:829
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_HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
virtual void v_GetEdgeToElementMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray)
Definition: StdTriExp.cpp:910
virtual const LibUtilities::BasisKey v_DetEdgeBasisKey(const int edge) const
Definition: StdTriExp.cpp:862
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space...