Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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  int P)
916  {
919  "Mapping not defined for this type of basis");
920 
921  int i;
922  int numModes;
923  int order0 = m_base[0]->GetNumModes();
924  int order1 = m_base[1]->GetNumModes();
925 
926  switch (eid)
927  {
928  case 0:
929  numModes = order0;
930  break;
931  case 1:
932  case 2:
933  numModes = order1;
934  break;
935  }
936 
937  bool checkForZeroedModes = false;
938  if (P == -1)
939  {
940  P = numModes;
941  }
942  else if(P != numModes)
943  {
944  checkForZeroedModes = true;
945  }
946 
947 
948  if(maparray.num_elements() != P)
949  {
950  maparray = Array<OneD, unsigned int>(P);
951  }
952 
953  if(signarray.num_elements() != P)
954  {
955  signarray = Array<OneD, int>(P,1);
956  }
957  else
958  {
959  fill(signarray.get() , signarray.get()+P, 1);
960  }
961 
962  switch(eid)
963  {
964  case 0:
965  {
966  int cnt = 0;
967  for(i = 0; i < P; cnt+=order1-i, ++i)
968  {
969  maparray[i] = cnt;
970  }
971 
972  if(edgeOrient==eBackwards)
973  {
974  swap( maparray[0] , maparray[1] );
975 
976  for(i = 3; i < P; i+=2)
977  {
978  signarray[i] = -1;
979  }
980  }
981  break;
982  }
983  case 1:
984  {
985  maparray[0] = order1;
986  maparray[1] = 1;
987  for(i = 2; i < P; i++)
988  {
989  maparray[i] = order1-1+i;
990  }
991 
992  if(edgeOrient==eBackwards)
993  {
994  swap( maparray[0] , maparray[1] );
995 
996  for(i = 3; i < P; i+=2)
997  {
998  signarray[i] = -1;
999  }
1000  }
1001  break;
1002  }
1003  case 2:
1004  {
1005  for(i = 0; i < P; i++)
1006  {
1007  maparray[i] = i;
1008  }
1009 
1010  if(edgeOrient==eForwards)
1011  {
1012  swap( maparray[0] , maparray[1] );
1013 
1014  for(i = 3; i < P; i+=2)
1015  {
1016  signarray[i] = -1;
1017  }
1018  }
1019  break;
1020  }
1021  default:
1022  ASSERTL0(false,"eid must be between 0 and 2");
1023  break;
1024  }
1025 
1026 
1027  if (checkForZeroedModes)
1028  {
1029  // Zero signmap and set maparray to zero if
1030  // elemental modes are not as large as face modes
1031  for (int j = numModes; j < P; j++)
1032  {
1033  signarray[j] = 0.0;
1034  maparray[j] = maparray[0];
1035  }
1036  }
1037  }
1038 
1039  int StdTriExp::v_GetVertexMap(const int localVertexId,bool useCoeffPacking)
1040  {
1041  ASSERTL0(
1042  GetEdgeBasisType(localVertexId) == LibUtilities::eModified_A ||
1043  GetEdgeBasisType(localVertexId) == LibUtilities::eModified_B,
1044  "Mapping not defined for this type of basis");
1045 
1046  int localDOF = 0;
1047  if(useCoeffPacking == true)
1048  {
1049  switch(localVertexId)
1050  {
1051  case 0:
1052  {
1053  localDOF = 0;
1054  break;
1055  }
1056  case 1:
1057  {
1058  localDOF = 1;
1059  break;
1060  }
1061  case 2:
1062  {
1063  localDOF = m_base[1]->GetNumModes();
1064  break;
1065  }
1066  default:
1067  {
1068  ASSERTL0(false,"eid must be between 0 and 2");
1069  break;
1070  }
1071  }
1072  }
1073  else // follow book format for vertex indexing.
1074  {
1075  switch(localVertexId)
1076  {
1077  case 0:
1078  {
1079  localDOF = 0;
1080  break;
1081  }
1082  case 1:
1083  {
1084  localDOF = m_base[1]->GetNumModes();
1085  break;
1086  }
1087  case 2:
1088  {
1089  localDOF = 1;
1090  break;
1091  }
1092  default:
1093  {
1094  ASSERTL0(false,"eid must be between 0 and 2");
1095  break;
1096  }
1097  }
1098  }
1099 
1100  return localDOF;
1101  }
1102 
1104  const int eid,
1105  const Orientation edgeOrient,
1106  Array<OneD, unsigned int>& maparray,
1107  Array<OneD, int>& signarray)
1108  {
1111  "Mapping not defined for this type of basis");
1112  int i;
1113  const int nummodes1 = m_base[1]->GetNumModes();
1114  const int nEdgeIntCoeffs = GetEdgeNcoeffs(eid)-2;
1115 
1116  if(maparray.num_elements() != nEdgeIntCoeffs)
1117  {
1118  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1119  }
1120 
1121  if(signarray.num_elements() != nEdgeIntCoeffs)
1122  {
1123  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
1124  }
1125  else
1126  {
1127  fill( signarray.get() , signarray.get()+nEdgeIntCoeffs, 1 );
1128  }
1129 
1130  switch(eid)
1131  {
1132  case 0:
1133  {
1134  int cnt = 2*nummodes1 - 1;
1135  for(i = 0; i < nEdgeIntCoeffs; cnt+=nummodes1-2-i, ++i)
1136  {
1137  maparray[i] = cnt;
1138  }
1139 
1140  if(edgeOrient==eBackwards)
1141  {
1142  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1143  {
1144  signarray[i] = -1;
1145  }
1146  }
1147  break;
1148  }
1149  case 1:
1150  {
1151  for(i = 0; i < nEdgeIntCoeffs; i++)
1152  {
1153  maparray[i] = nummodes1+1+i;
1154  }
1155 
1156  if(edgeOrient==eBackwards)
1157  {
1158  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1159  {
1160  signarray[i] = -1;
1161  }
1162  }
1163  break;
1164  }
1165  case 2:
1166  {
1167  for(i = 0; i < nEdgeIntCoeffs; i++)
1168  {
1169  maparray[i] = 2+i;
1170  }
1171 
1172  if(edgeOrient==eForwards)
1173  {
1174  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1175  {
1176  signarray[i] = -1;
1177  }
1178  }
1179  break;
1180  }
1181  default:
1182  {
1183  ASSERTL0(false,"eid must be between 0 and 2");
1184  break;
1185  }
1186  }
1187  }
1188 
1190  {
1193  "Expansion not of a proper type");
1194 
1195  int i,j;
1196  int cnt = 0;
1197  int nummodes0, nummodes1;
1198  int startvalue;
1199  if(outarray.num_elements()!=GetNcoeffs()-NumBndryCoeffs())
1200  {
1202  }
1203 
1204  nummodes0 = m_base[0]->GetNumModes();
1205  nummodes1 = m_base[1]->GetNumModes();
1206 
1207  startvalue = 2*nummodes1;
1208 
1209  for(i = 0; i < nummodes0-2; i++)
1210  {
1211  for(j = 0; j < nummodes1-3-i; j++)
1212  {
1213  outarray[cnt++]=startvalue+j;
1214  }
1215  startvalue+=nummodes1-2-i;
1216  }
1217  }
1218 
1220  {
1223  "Expansion not of expected type");
1224  int i;
1225  int cnt;
1226  int nummodes0, nummodes1;
1227  int value;
1228 
1229  if (outarray.num_elements()!=NumBndryCoeffs())
1230  {
1232  }
1233 
1234  nummodes0 = m_base[0]->GetNumModes();
1235  nummodes1 = m_base[1]->GetNumModes();
1236 
1237  value = 2*nummodes1-1;
1238  for(i = 0; i < value; i++)
1239  {
1240  outarray[i]=i;
1241  }
1242  cnt = value;
1243 
1244  for(i = 0; i < nummodes0-2; i++)
1245  {
1246  outarray[cnt++]=value;
1247  value += nummodes1-2-i;
1248  }
1249  }
1250 
1251 
1252  //---------------------------------------
1253  // Wrapper functions
1254  //---------------------------------------
1255 
1257  {
1258 
1259  MatrixType mtype = mkey.GetMatrixType();
1260 
1261  DNekMatSharedPtr Mat;
1262 
1263  switch(mtype)
1264  {
1266  {
1267  int nq0, nq1, nq;
1268 
1269  nq0 = m_base[0]->GetNumPoints();
1270  nq1 = m_base[1]->GetNumPoints();
1271 
1272  // take definition from key
1274  {
1275  nq = (int) mkey.GetConstFactor(eFactorConst);
1276  }
1277  else
1278  {
1279  nq = max(nq0,nq1);
1280  }
1281 
1282  int neq = LibUtilities::StdTriData::
1284  Array<OneD, Array<OneD, NekDouble> > coords(neq);
1285  Array<OneD, NekDouble> coll (2);
1287  Array<OneD, NekDouble> tmp (nq0);
1288 
1289 
1290  Mat = MemoryManager<DNekMat>::AllocateSharedPtr(neq,nq0*nq1);
1291  int cnt = 0;
1292 
1293  for(int i = 0; i < nq; ++i)
1294  {
1295  for(int j = 0; j < nq-i; ++j,++cnt)
1296  {
1297  coords[cnt] = Array<OneD, NekDouble>(2);
1298  coords[cnt][0] = -1.0 + 2*j/(NekDouble)(nq-1);
1299  coords[cnt][1] = -1.0 + 2*i/(NekDouble)(nq-1);
1300  }
1301  }
1302 
1303  for(int i = 0; i < neq; ++i)
1304  {
1305  LocCoordToLocCollapsed(coords[i],coll);
1306 
1307  I[0] = m_base[0]->GetI(coll);
1308  I[1] = m_base[1]->GetI(coll+1);
1309 
1310  // interpolate first coordinate direction
1311  for (int j = 0; j < nq1; ++j)
1312  {
1313  NekDouble fac = (I[1]->GetPtr())[j];
1314  Vmath::Smul(nq0, fac, I[0]->GetPtr(), 1, tmp, 1);
1315 
1316  Vmath::Vcopy(nq0, &tmp[0], 1,
1317  Mat->GetRawPtr() + j*nq0*neq + i, neq);
1318  }
1319 
1320  }
1321  break;
1322  }
1323  default:
1324  {
1326  break;
1327  }
1328  }
1329 
1330  return Mat;
1331  }
1332 
1334  {
1335  return v_GenMatrix(mkey);
1336  }
1337 
1338 
1339  //---------------------------------------
1340  // Operator evaluation functions
1341  //---------------------------------------
1342 
1344  const Array<OneD, const NekDouble> &inarray,
1345  Array<OneD, NekDouble> &outarray,
1346  const StdMatrixKey &mkey)
1347  {
1348  StdExpansion::MassMatrixOp_MatFree(inarray,outarray,mkey);
1349  }
1350 
1352  const Array<OneD, const NekDouble> &inarray,
1353  Array<OneD, NekDouble> &outarray,
1354  const StdMatrixKey &mkey)
1355  {
1356  StdTriExp::v_LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
1357  }
1358 
1360  const int k1,
1361  const int k2,
1362  const Array<OneD, const NekDouble> &inarray,
1363  Array<OneD, NekDouble> &outarray,
1364  const StdMatrixKey &mkey)
1365  {
1367  k1,k2,inarray,outarray,mkey);
1368  }
1369 
1371  const int i,
1372  const Array<OneD, const NekDouble> &inarray,
1373  Array<OneD, NekDouble> &outarray,
1374  const StdMatrixKey &mkey)
1375  {
1376  StdExpansion::WeakDerivMatrixOp_MatFree(i,inarray,outarray,mkey);
1377  }
1378 
1380  const Array<OneD, const NekDouble> &inarray,
1381  Array<OneD, NekDouble> &outarray,
1382  const StdMatrixKey &mkey)
1383  {
1384  StdTriExp::v_HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
1385  }
1386 
1387 
1389  const StdMatrixKey &mkey)
1390  {
1391  int qa = m_base[0]->GetNumPoints();
1392  int qb = m_base[1]->GetNumPoints();
1393  int nmodes_a = m_base[0]->GetNumModes();
1394  int nmodes_b = m_base[1]->GetNumModes();
1395  int nmodes = min(nmodes_a,nmodes_b);
1396 
1397  // Declare orthogonal basis.
1400 
1403  StdTriExp OrthoExp(Ba,Bb);
1404 
1405  NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDiffCoeff);
1406 
1407  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
1408 
1409 
1410  if(mkey.HasVarCoeff(eVarCoeffLaplacian)) // Rodrigo's svv mapping
1411  {
1412  Array<OneD, NekDouble> sqrt_varcoeff(qa*qb);
1413  Array<OneD, NekDouble> tmp(qa*qb);
1414 
1415  Vmath::Vsqrt(qa * qb,
1417  sqrt_varcoeff, 1);
1418 
1419  // multiply by sqrt(Variable Coefficient) containing h v /p
1420  Vmath::Vmul(qa*qb,sqrt_varcoeff,1,array,1,tmp,1);
1421 
1422  // project onto modal space.
1423  OrthoExp.FwdTrans(tmp,orthocoeffs);
1424 
1425  int cnt = 0;
1426  for(int j = 0; j < nmodes_a; ++j)
1427  {
1428  for(int k = 0; k < nmodes_b-j; ++k, ++cnt)
1429  {
1430  orthocoeffs[cnt] *=
1431  (1.0 + SvvDiffCoeff
1432  *pow(j/(nmodes_a-1)+k/(nmodes_b-1),0.5*nmodes));
1433  }
1434  }
1435 
1436  // backward transform to physical space
1437  OrthoExp.BwdTrans(orthocoeffs,tmp);
1438 
1439  // multiply by sqrt(Variable Coefficient) containing h v /p
1440  // - split to keep symmetry
1441  Vmath::Vmul(qa*qb,sqrt_varcoeff,1,tmp,1,array,1);
1442  }
1443  else
1444  {
1445  int j, k , cnt = 0;
1446  int cutoff = (int) (mkey.GetConstFactor(eFactorSVVCutoffRatio)*
1447  min(nmodes_a,nmodes_b));
1448 
1449  NekDouble epsilon = 1.0;
1450  int nmodes = min(nmodes_a,nmodes_b);
1451 
1452  // project onto physical space.
1453  OrthoExp.FwdTrans(array,orthocoeffs);
1454 
1455  // apply SVV filter (JEL)
1456  for(j = 0; j < nmodes_a; ++j)
1457  {
1458  for(k = 0; k < nmodes_b-j; ++k)
1459  {
1460  if(j + k >= cutoff)
1461  {
1462  orthocoeffs[cnt] *= (SvvDiffCoeff
1463  *exp(-(j+k-nmodes)*(j+k-nmodes)
1464  /((NekDouble)((j+k-cutoff+epsilon)
1465  *(j+k-cutoff+epsilon)))));
1466  }
1467  else
1468  {
1469  orthocoeffs[cnt] *= 0.0;
1470  }
1471  cnt++;
1472  }
1473  }
1474 
1475  // backward transform to physical space
1476  OrthoExp.BwdTrans(orthocoeffs,array);
1477  }
1478  }
1479 
1481  int numMin,
1482  const Array<OneD, const NekDouble> &inarray,
1483  Array<OneD, NekDouble> &outarray)
1484  {
1485  int n_coeffs = inarray.num_elements();
1486  int nquad0 = m_base[0]->GetNumPoints();
1487  int nquad1 = m_base[1]->GetNumPoints();
1488  Array<OneD, NekDouble> coeff(n_coeffs);
1489  Array<OneD, NekDouble> coeff_tmp(n_coeffs,0.0);
1492  int nqtot = nquad0*nquad1;
1493  Array<OneD, NekDouble> phys_tmp(nqtot,0.0);
1494 
1495  int nmodes0 = m_base[0]->GetNumModes();
1496  int nmodes1 = m_base[1]->GetNumModes();
1497  int numMin2 = nmodes0;
1498  int i;
1499 
1500  const LibUtilities::PointsKey Pkey0(
1502  const LibUtilities::PointsKey Pkey1(
1504 
1505  LibUtilities::BasisKey b0(m_base[0]->GetBasisType(),nmodes0,Pkey0);
1506  LibUtilities::BasisKey b1(m_base[1]->GetBasisType(),nmodes1,Pkey1);
1507 
1508  LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A,nmodes0,Pkey0);
1509  LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_B,nmodes1,Pkey1);
1510 
1511  StdRegions::StdTriExpSharedPtr m_OrthoTriExp;
1513 
1515  ::AllocateSharedPtr(b0, b1);
1516  m_OrthoTriExp = MemoryManager<StdRegions::StdTriExp>
1517  ::AllocateSharedPtr(bortho0, bortho1);
1518 
1519  m_TriExp ->BwdTrans(inarray,phys_tmp);
1520  m_OrthoTriExp->FwdTrans(phys_tmp, coeff);
1521 
1522  for (i = 0; i < n_coeffs; i++)
1523  {
1524  if (i == numMin)
1525  {
1526  coeff[i] = 0.0;
1527  numMin += numMin2 - 1;
1528  numMin2 -= 1.0;
1529  }
1530  }
1531 
1532  m_OrthoTriExp->BwdTrans(coeff,phys_tmp);
1533  m_TriExp ->FwdTrans(phys_tmp, outarray);
1534 
1535  }
1536 
1538  const Array<OneD, const NekDouble> &inarray,
1539  Array<OneD, NekDouble> &outarray,
1540  const StdMatrixKey &mkey)
1541  {
1543 
1544  if(inarray.get() == outarray.get())
1545  {
1547  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
1548 
1549  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, 1.0, mat->GetPtr().get(),
1550  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
1551  }
1552  else
1553  {
1554  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, 1.0, mat->GetPtr().get(),
1555  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
1556  }
1557  }
1558 
1559  //---------------------------------------
1560  // Private helper functions
1561  //---------------------------------------
1562 
1564  const Array<OneD, const NekDouble>& inarray,
1565  Array<OneD, NekDouble> &outarray)
1566  {
1567  int i;
1568  int nquad0 = m_base[0]->GetNumPoints();
1569  int nquad1 = m_base[1]->GetNumPoints();
1570 
1571  const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
1572  const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
1573  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
1574 
1575  // multiply by integration constants
1576  for(i = 0; i < nquad1; ++i)
1577  {
1578  Vmath::Vmul(nquad0,inarray.get()+i*nquad0,1,
1579  w0.get(),1, outarray.get()+i*nquad0,1);
1580  }
1581 
1582  switch(m_base[1]->GetPointsType())
1583  {
1584  // Legendre inner product
1587  for(i = 0; i < nquad1; ++i)
1588  {
1589  Blas::Dscal(nquad0,0.5*(1-z1[i])*w1[i],
1590  outarray.get()+i*nquad0,1);
1591  }
1592  break;
1593 
1594  // (1,0) Jacobi Inner product
1596  for(i = 0; i < nquad1; ++i)
1597  {
1598  Blas::Dscal(nquad0,0.5*w1[i],outarray.get()+i*nquad0,1);
1599  }
1600  break;
1601 
1602  default:
1603  ASSERTL0(false, "Unsupported quadrature points type.");
1604  break;
1605  }
1606  }
1607 
1609  Array<OneD, int> &conn,
1610  bool standard)
1611  {
1612  int np1 = m_base[0]->GetNumPoints();
1613  int np2 = m_base[1]->GetNumPoints();
1614  int np = max(np1,np2);
1615 
1616  conn = Array<OneD, int>(3*(np-1)*(np-1));
1617 
1618  int row = 0;
1619  int rowp1 = 0;
1620  int cnt = 0;
1621  for(int i = 0; i < np-1; ++i)
1622  {
1623  rowp1 += np-i;
1624  for(int j = 0; j < np-i-2; ++j)
1625  {
1626  conn[cnt++] = row +j;
1627  conn[cnt++] = row +j+1;
1628  conn[cnt++] = rowp1 +j;
1629 
1630  conn[cnt++] = rowp1 +j+1;
1631  conn[cnt++] = rowp1 +j;
1632  conn[cnt++] = row +j+1;
1633  }
1634 
1635  conn[cnt++] = row +np-i-2;
1636  conn[cnt++] = row +np-i-1;
1637  conn[cnt++] = rowp1+np-i-2;
1638 
1639  row += np-i;
1640  }
1641  }
1642 
1643 
1644  }//end namespace
1645 }//end namespace
boost::shared_ptr< StdTriExp > StdTriExpSharedPtr
Definition: StdTriExp.h:267
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:470
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:1537
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
virtual void v_GetEdgeToElementMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, int P=-1)
Definition: StdTriExp.cpp:910
void MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:971
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
bool HasVarCoeff(const StdRegions::VarCoeffType &coeff) const
Definition: StdMatrixKey.h:173
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:22
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:394
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:942
virtual void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:1563
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:1370
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
const Array< OneD, const NekDouble > & GetVarCoeff(const StdRegions::VarCoeffType &coeff) const
Definition: StdMatrixKey.h:152
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:1379
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:1189
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
Definition: StdTriExp.cpp:1388
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:700
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:1256
Principle Orthogonal Functions .
Definition: BasisType.h:47
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:131
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:1608
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:1480
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:413
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:1219
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:1343
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:1351
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
Definition: StdTriExp.cpp:1333
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:525
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
Definition: StdTriExp.cpp:1039
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:1103
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:1047
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 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...