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