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