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