Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
StdQuadExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File StdQuadExp.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: Quadrilateral routines built upon StdExpansion2D
33 //
34 ///////////////////////////////////////////////////////////////////////////////
36 #include <StdRegions/StdQuadExp.h>
37 #include <StdRegions/StdSegExp.h>
39 
40 using namespace std;
41 
42 namespace Nektar
43 {
44  namespace StdRegions
45  {
46 
47 
48  StdQuadExp::StdQuadExp()
49  {
50  }
51 
52  /** \brief Constructor using BasisKey class for quadrature
53  * points and order definition
54  */
55  StdQuadExp::StdQuadExp(const LibUtilities::BasisKey &Ba,
56  const LibUtilities::BasisKey &Bb):
57  StdExpansion (Ba.GetNumModes()*Bb.GetNumModes(),2,Ba,Bb),
58  StdExpansion2D(Ba.GetNumModes()*Bb.GetNumModes(),Ba,Bb)
59  {
60  }
61 
62  /** \brief Copy Constructor */
64  StdExpansion(T),
66  {
67  }
68 
69  /** \brief Destructor */
71  {
72  }
73 
74  /////////////////////////
75  // Integration Methods //
76  /////////////////////////
77 
79  const Array<OneD, const NekDouble>& inarray)
80  {
81  Array<OneD, const NekDouble> w0 = m_base[0]->GetW();
82  Array<OneD, const NekDouble> w1 = m_base[1]->GetW();
83 
84  return StdExpansion2D::Integral(inarray,w0,w1);
85  }
86 
87  /////////////////////////////
88  // Differentiation Methods //
89  /////////////////////////////
90 
91  /** \brief Calculate the derivative of the physical points
92  *
93  * For quadrilateral region can use the Tensor_Deriv function
94  * defined under StdExpansion.
95  */
96 
98  Array<OneD, NekDouble> &out_d0,
99  Array<OneD, NekDouble> &out_d1,
100  Array<OneD, NekDouble> &out_d2)
101  {
102  PhysTensorDeriv(inarray, out_d0, out_d1);
103  }
104 
105  void StdQuadExp::v_PhysDeriv(const int dir,
106  const Array<OneD, const NekDouble>& inarray,
107  Array<OneD, NekDouble> &outarray)
108  {
109  switch(dir)
110  {
111  case 0:
112  {
113  PhysTensorDeriv(inarray, outarray, NullNekDouble1DArray);
114  }
115  break;
116  case 1:
117  {
118  PhysTensorDeriv(inarray, NullNekDouble1DArray, outarray);
119  }
120  break;
121  default:
122  {
123  ASSERTL1(false,"input dir is out of range");
124  }
125  break;
126  }
127  }
128 
130  Array<OneD, NekDouble> &out_d0,
131  Array<OneD, NekDouble> &out_d1,
132  Array<OneD, NekDouble> &out_d2)
133  {
134  //PhysTensorDeriv(inarray, out_d0, out_d1);
135  StdQuadExp::v_PhysDeriv(inarray, out_d0, out_d1);
136  }
137 
138  void StdQuadExp::v_StdPhysDeriv(const int dir,
139  const Array<OneD, const NekDouble>& inarray,
140  Array<OneD, NekDouble> &outarray)
141  {
142  //PhysTensorDeriv(inarray, outarray);
143  StdQuadExp::v_PhysDeriv(dir,inarray,outarray);
144 
145  }
146 
147 
148  ////////////////
149  // Transforms //
150  ////////////////
151 
153  Array<OneD, NekDouble> &outarray)
154  {
155  if(m_base[0]->Collocation() && m_base[1]->Collocation())
156  {
158  inarray, 1, outarray, 1);
159  }
160  else
161  {
162  StdQuadExp::v_BwdTrans_SumFac(inarray,outarray);
163  }
164  }
165 
167  const Array<OneD, const NekDouble>& inarray,
168  Array<OneD, NekDouble> &outarray)
169  {
171  m_base[1]->GetNumModes());
172 
173  BwdTrans_SumFacKernel(m_base[0]->GetBdata(),
174  m_base[1]->GetBdata(),
175  inarray,outarray,wsp,true,true);
176  }
177 
178  // The arguments doCheckCollDir0 and doCheckCollDir1 allow you to specify whether
179  // to check if the basis has collocation properties (i.e. for the classical spectral
180  // element basis, In this case the 1D 'B' matrix is equal to the identity matrix
181  // which can be exploited to speed up the calculations).
182  // However, as this routine also allows to pass the matrix 'DB' (derivative of the basis),
183  // the collocation property cannot always be used. Therefor follow this rule:
184  // if base0 == m_base[0]->GetBdata() --> set doCheckCollDir0 == true;
185  // base1 == m_base[1]->GetBdata() --> set doCheckCollDir1 == true;
186  // base0 == m_base[0]->GetDbdata() --> set doCheckCollDir0 == false;
187  // base1 == m_base[1]->GetDbdata() --> set doCheckCollDir1 == false;
189  const Array<OneD, const NekDouble>& base0,
190  const Array<OneD, const NekDouble>& base1,
191  const Array<OneD, const NekDouble>& inarray,
192  Array<OneD, NekDouble> &outarray,
194  bool doCheckCollDir0,
195  bool doCheckCollDir1)
196  {
197  int nquad0 = m_base[0]->GetNumPoints();
198  int nquad1 = m_base[1]->GetNumPoints();
199  int nmodes0 = m_base[0]->GetNumModes();
200  int nmodes1 = m_base[1]->GetNumModes();
201 
202  bool colldir0 = doCheckCollDir0?(m_base[0]->Collocation()):false;
203  bool colldir1 = doCheckCollDir1?(m_base[1]->Collocation()):false;
204 
205  if(colldir0 && colldir1)
206  {
207  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,outarray.get(),1);
208  }
209  else if(colldir0)
210  {
211  Blas::Dgemm('N','T', nquad0, nquad1,nmodes1, 1.0, &inarray[0], nquad0,
212  base1.get(), nquad1, 0.0, &outarray[0], nquad0);
213  }
214  else if(colldir1)
215  {
216  Blas::Dgemm('N','N', nquad0,nmodes1,nmodes0,1.0, base0.get(),
217  nquad0, &inarray[0], nmodes0,0.0,&outarray[0], nquad0);
218  }
219  else
220  {
221  ASSERTL1(wsp.num_elements()>=nquad0*nmodes1,"Workspace size is not sufficient");
222 
223  // Those two calls correpsond to the operation
224  // out = B0*in*Transpose(B1);
225  Blas::Dgemm('N','N', nquad0,nmodes1,nmodes0,1.0, base0.get(),
226  nquad0, &inarray[0], nmodes0,0.0,&wsp[0], nquad0);
227  Blas::Dgemm('N','T', nquad0, nquad1,nmodes1, 1.0, &wsp[0], nquad0,
228  base1.get(), nquad1, 0.0, &outarray[0], nquad0);
229  }
230  }
231 
232 
234  Array<OneD, NekDouble> &outarray)
235  {
236  if((m_base[0]->Collocation())&&(m_base[1]->Collocation()))
237  {
238  Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
239  }
240  else
241  {
242  StdQuadExp::v_IProductWRTBase(inarray,outarray);
243 
244  // get Mass matrix inverse
245  StdMatrixKey masskey(eInvMass,DetShapeType(),*this);
246  DNekMatSharedPtr matsys = GetStdMatrix(masskey);
247 
248  // copy inarray in case inarray == outarray
249  NekVector<NekDouble> in(m_ncoeffs,outarray,eCopy);
251 
252  out = (*matsys)*in;
253  }
254  }
255 
257  const Array<OneD, const NekDouble>& inarray,
258  Array<OneD, NekDouble> &outarray)
259  {
260  if((m_base[0]->Collocation())&&(m_base[1]->Collocation()))
261  {
262  Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
263  }
264  else
265  {
266  int i,j;
267  int npoints[2] = {m_base[0]->GetNumPoints(),
268  m_base[1]->GetNumPoints()};
269  int nmodes[2] = {m_base[0]->GetNumModes(),
270  m_base[1]->GetNumModes()};
271 
272  fill(outarray.get(), outarray.get()+m_ncoeffs, 0.0 );
273 
274  Array<OneD, NekDouble> physEdge[4];
275  Array<OneD, NekDouble> coeffEdge[4];
276  for(i = 0; i < 4; i++)
277  {
278  physEdge[i] = Array<OneD, NekDouble>(npoints[i%2]);
279  coeffEdge[i] = Array<OneD, NekDouble>(nmodes[i%2]);
280  }
281 
282  for(i = 0; i < npoints[0]; i++)
283  {
284  physEdge[0][i] = inarray[i];
285  physEdge[2][i] = inarray[npoints[0]*npoints[1]-1-i];
286  }
287 
288  for(i = 0; i < npoints[1]; i++)
289  {
290  physEdge[1][i] = inarray[npoints[0]-1+i*npoints[0]];
291  physEdge[3][i] = inarray[(npoints[1]-1)*npoints[0]-i*npoints[0]];
292  }
293 
296 
297  Array<OneD, unsigned int> mapArray;
298  Array<OneD, int> signArray;
299  NekDouble sign;
300 
301  for(i = 0; i < 4; i++)
302  {
303  segexp[i%2]->FwdTrans_BndConstrained(physEdge[i],coeffEdge[i]);
304 
305  GetEdgeToElementMap(i,eForwards,mapArray,signArray);
306  for(j=0; j < nmodes[i%2]; j++)
307  {
308  sign = (NekDouble) signArray[j];
309  outarray[ mapArray[j] ] = sign * coeffEdge[i][j];
310  }
311  }
312 
315 
316  StdMatrixKey masskey(eMass,DetShapeType(),*this);
317  MassMatrixOp(outarray,tmp0,masskey);
318  IProductWRTBase(inarray,tmp1);
319 
320  Vmath::Vsub(m_ncoeffs, tmp1, 1, tmp0, 1, tmp1, 1);
321 
322  // get Mass matrix inverse (only of interior DOF)
323  // use block (1,1) of the static condensed system
324  // note: this block alreay contains the inverse matrix
325  DNekMatSharedPtr matsys = (m_stdStaticCondMatrixManager[masskey])->GetBlock(1,1);
326 
327  int nBoundaryDofs = NumBndryCoeffs();
328  int nInteriorDofs = m_ncoeffs - nBoundaryDofs;
329 
330 
331  Array<OneD, NekDouble> rhs(nInteriorDofs);
332  Array<OneD, NekDouble> result(nInteriorDofs);
333 
334  GetInteriorMap(mapArray);
335 
336  for(i = 0; i < nInteriorDofs; i++)
337  {
338  rhs[i] = tmp1[ mapArray[i] ];
339  }
340 
341  Blas::Dgemv('N',nInteriorDofs,nInteriorDofs,1.0, &(matsys->GetPtr())[0],
342  nInteriorDofs,rhs.get(),1,0.0,result.get(),1);
343 
344  for(i = 0; i < nInteriorDofs; i++)
345  {
346  outarray[ mapArray[i] ] = result[i];
347  }
348  }
349 
350  }
351 
352  /////////////////////////////
353  // Inner Product Functions //
354  /////////////////////////////
355 
356  /** \brief Calculate the inner product of inarray with respect to
357  * the basis B=base0*base1 and put into outarray
358  *
359  * \f$
360  * \begin{array}{rcl}
361  * I_{pq} = (\phi_q \phi_q, u) & = & \sum_{i=0}^{nq_0}
362  * \sum_{j=0}^{nq_1}
363  * \phi_p(\xi_{0,i}) \phi_q(\xi_{1,j}) w^0_i w^1_j u(\xi_{0,i}
364  * \xi_{1,j}) \\
365  * & = & \sum_{i=0}^{nq_0} \phi_p(\xi_{0,i})
366  * \sum_{j=0}^{nq_1} \phi_q(\xi_{1,j}) \tilde{u}_{i,j}
367  * \end{array}
368  * \f$
369  *
370  * where
371  *
372  * \f$ \tilde{u}_{i,j} = w^0_i w^1_j u(\xi_{0,i},\xi_{1,j}) \f$
373  *
374  * which can be implemented as
375  *
376  * \f$ f_{qi} = \sum_{j=0}^{nq_1} \phi_q(\xi_{1,j})
377  * \tilde{u}_{i,j} = {\bf B_1 U} \f$
378  * \f$ I_{pq} = \sum_{i=0}^{nq_0} \phi_p(\xi_{0,i}) f_{qi} =
379  * {\bf B_0 F} \f$
380  */
382  const Array<OneD, const NekDouble>& inarray,
383  Array<OneD, NekDouble> &outarray)
384  {
385  if(m_base[0]->Collocation() && m_base[1]->Collocation())
386  {
387  MultiplyByQuadratureMetric(inarray,outarray);
388  }
389  else
390  {
391  StdQuadExp::v_IProductWRTBase_SumFac(inarray,outarray);
392  }
393  }
394 
396  const Array<OneD, const NekDouble>& inarray,
397  Array<OneD, NekDouble> &outarray,
398  bool multiplybyweights)
399  {
400  int nquad0 = m_base[0]->GetNumPoints();
401  int nquad1 = m_base[1]->GetNumPoints();
402  int order0 = m_base[0]->GetNumModes();
403 
404  if(multiplybyweights)
405  {
406  Array<OneD,NekDouble> tmp(nquad0*nquad1+nquad1*order0);
407  Array<OneD,NekDouble> wsp(tmp+nquad0*nquad1);
408 
409  // multiply by integration constants
410  MultiplyByQuadratureMetric(inarray,tmp);
411  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
412  m_base[1]->GetBdata(),
413  tmp,outarray,wsp,true,true);
414  }
415  else
416  {
417  Array<OneD,NekDouble> wsp(nquad1*order0);
418  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
419  m_base[1]->GetBdata(),
420  inarray,outarray,wsp,true,true);
421  }
422  }
423 
425  const Array<OneD, const NekDouble>& inarray,
426  Array<OneD, NekDouble> &outarray)
427  {
428  int nq = GetTotPoints();
429  StdMatrixKey iprodmatkey(eIProductWRTBase,DetShapeType(),*this);
430  DNekMatSharedPtr iprodmat = GetStdMatrix(iprodmatkey);
431 
432  Blas::Dgemv('N',m_ncoeffs,nq,1.0,iprodmat->GetPtr().get(),
433  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
434  }
435 
437  const Array<OneD, const NekDouble>& inarray,
438  Array<OneD, NekDouble> & outarray)
439  {
440  StdQuadExp::IProductWRTDerivBase_SumFac(dir,inarray,outarray);
441  }
442 
444  const Array<OneD, const NekDouble>& inarray,
445  Array<OneD, NekDouble> &outarray)
446  {
447  ASSERTL0((dir==0)||(dir==1),"input dir is out of range");
448 
449  int nquad0 = m_base[0]->GetNumPoints();
450  int nquad1 = m_base[1]->GetNumPoints();
451  int nqtot = nquad0*nquad1;
452  int order0 = m_base[0]->GetNumModes();
453 
454  Array<OneD,NekDouble> tmp(nqtot+nquad1*order0);
455  Array<OneD,NekDouble> wsp(tmp+nqtot);
456 
457  // multiply by integration constants
458  MultiplyByQuadratureMetric(inarray,tmp);
459 
460  if(dir) // dir == 1
461  {
462  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
463  m_base[1]->GetDbdata(),
464  tmp,outarray,wsp,true,false);
465  }
466  else // dir == 0
467  {
468  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
469  m_base[1]->GetBdata(),
470  tmp,outarray,wsp,false,true);
471  }
472  }
473 
475  const Array<OneD, const NekDouble>& inarray,
476  Array<OneD, NekDouble> &outarray)
477  {
478  ASSERTL0((dir==0)||(dir==1),"input dir is out of range");
479 
480  int nq = GetTotPoints();
481  MatrixType mtype;
482 
483  if(dir) // dir == 1
484  {
485  mtype = eIProductWRTDerivBase1;
486  }
487  else // dir == 0
488  {
489  mtype = eIProductWRTDerivBase0;
490  }
491 
492  StdMatrixKey iprodmatkey(mtype,DetShapeType(),*this);
493  DNekMatSharedPtr iprodmat = GetStdMatrix(iprodmatkey);
494 
495  Blas::Dgemv('N',m_ncoeffs,nq,1.0,iprodmat->GetPtr().get(),
496  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
497  }
498 
499  // the arguments doCheckCollDir0 and doCheckCollDir1 allow you to specify whether
500  // to check if the basis has collocation properties (i.e. for the classical spectral
501  // element basis, In this case the 1D 'B' matrix is equal to the identity matrix
502  // which can be exploited to speed up the calculations).
503  // However, as this routine also allows to pass the matrix 'DB' (derivative of the basis),
504  // the collocation property cannot always be used. Therefor follow this rule:
505  // if base0 == m_base[0]->GetBdata() --> set doCheckCollDir0 == true;
506  // base1 == m_base[1]->GetBdata() --> set doCheckCollDir1 == true;
507  // base0 == m_base[0]->GetDbdata() --> set doCheckCollDir0 == false;
508  // base1 == m_base[1]->GetDbdata() --> set doCheckCollDir1 == false;
510  const Array<OneD, const NekDouble>& base0,
511  const Array<OneD, const NekDouble>& base1,
512  const Array<OneD, const NekDouble>& inarray,
513  Array<OneD, NekDouble> &outarray,
515  bool doCheckCollDir0,
516  bool doCheckCollDir1)
517  {
518  int nquad0 = m_base[0]->GetNumPoints();
519  int nquad1 = m_base[1]->GetNumPoints();
520  int nmodes0 = m_base[0]->GetNumModes();
521  int nmodes1 = m_base[1]->GetNumModes();
522 
523  bool colldir0 = doCheckCollDir0?(m_base[0]->Collocation()):false;
524  bool colldir1 = doCheckCollDir1?(m_base[1]->Collocation()):false;
525 
526  if(colldir0 && colldir1)
527  {
528  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,outarray.get(),1);
529  }
530  else if(colldir0)
531  {
532  Blas::Dgemm('N', 'N',nmodes0,nmodes1, nquad1,1.0, inarray.get(),
533  nmodes0, base1.get(), nquad1, 0.0,outarray.get(),nmodes0);
534  }
535  else if(colldir1)
536  {
537  Blas::Dgemm('T','N',nmodes0,nquad1,nquad0,1.0,base0.get(),
538  nquad0,inarray.get(),nquad0,0.0,outarray.get(),nmodes0);
539  }
540  else
541  {
542  ASSERTL1(wsp.num_elements()>=nquad1*nmodes0,"Workspace size is not sufficient");
543 
544 #if 1
545  Blas::Dgemm('T','N',nmodes0,nquad1,nquad0,1.0,base0.get(),
546  nquad0,inarray.get(),nquad0,0.0,wsp.get(),nmodes0);
547 
548 #else
549  for(int i = 0; i < nmodes0; ++i)
550  {
551  for(int j = 0; j < nquad1; ++j)
552  {
553  wsp[j*nmodes0+i] = Blas::Ddot(nquad0,
554  base0.get()+i*nquad0,1,
555  inarray.get()+j*nquad0,1);
556  }
557  }
558 #endif
559  Blas::Dgemm('N', 'N',nmodes0,nmodes1, nquad1,1.0, wsp.get(),
560  nmodes0, base1.get(), nquad1, 0.0,outarray.get(),nmodes0);
561  }
562  }
563 
564  //////////////////////////
565  // Evaluation functions //
566  //////////////////////////
567 
568 
571  {
572  eta[0] = xi[0];
573  eta[1] = xi[1];
574  }
575 
576  /** \brief Fill outarray with mode \a mode of expansion
577  *
578  * Note for quadrilateral expansions _base[0] (i.e. p) modes run
579  * fastest
580  */
581 
582  void StdQuadExp::v_FillMode(const int mode,
583  Array<OneD, NekDouble> &outarray)
584  {
585  int i;
586  int nquad0 = m_base[0]->GetNumPoints();
587  int nquad1 = m_base[1]->GetNumPoints();
588  Array<OneD, const NekDouble> base0 = m_base[0]->GetBdata();
589  Array<OneD, const NekDouble> base1 = m_base[1]->GetBdata();
590  int btmp0 = m_base[0]->GetNumModes();
591  int mode0 = mode%btmp0;
592  int mode1 = mode/btmp0;
593 
594 
595  ASSERTL2(mode1 == (int)floor((1.0*mode)/btmp0),
596  "Integer Truncation not Equiv to Floor");
597 
598  ASSERTL2(m_ncoeffs <= mode,
599  "calling argument mode is larger than total expansion order");
600 
601  for(i = 0; i < nquad1; ++i)
602  {
603  Vmath::Vcopy(nquad0,(NekDouble *)(base0.get() + mode0*nquad0),
604  1, &outarray[0]+i*nquad0,1);
605  }
606 
607  for(i = 0; i < nquad0; ++i)
608  {
609  Vmath::Vmul(nquad1,(NekDouble *)(base1.get() + mode1*nquad1),1,
610  &outarray[0]+i,nquad0,&outarray[0]+i,nquad0);
611  }
612  }
613 
614  //////////////////////
615  // Helper functions //
616  //////////////////////
617 
619  {
620  return 4;
621  }
622 
624  {
625  return 4;
626  }
627 
628  int StdQuadExp::v_GetEdgeNcoeffs(const int i) const
629  {
630  ASSERTL2((i >= 0)&&(i <= 3),"edge id is out of range");
631 
632  if((i == 0)||(i == 2))
633  {
634  return GetBasisNumModes(0);
635  }
636  else
637  {
638  return GetBasisNumModes(1);
639  }
640  }
641 
642  int StdQuadExp::v_GetEdgeNumPoints(const int i) const
643  {
644  ASSERTL2((i >= 0)&&(i <= 3),"edge id is out of range");
645 
646  if((i == 0)||(i == 2))
647  {
648  return GetNumPoints(0);
649  }
650  else
651  {
652  return GetNumPoints(1);
653  }
654  }
655 
657  {
658  ASSERTL2((i >= 0)&&(i <= 3),"edge id is out of range");
659 
660  if((i == 0)||(i == 2))
661  {
662  return GetBasisType(0);
663  }
664  else
665  {
666  return GetBasisType(1);
667  }
668  }
669 
671  {
672  ASSERTL2((edge >= 0)&&(edge <= 3),"edge id is out of range");
673 
674  if((edge == 0)||(edge == 2))
675  {
676  return 0;
677  }
678  else
679  {
680  return 1;
681  }
682  }
683 
685  {
686  ASSERTL2((i >= 0)&&(i <= 3),"edge id is out of range");
687 
688  if((i == 0)||(i == 2))
689  {
690  return GetBasis(0)->GetBasisKey();
691  }
692  else
693  {
694  return GetBasis(1)->GetBasisKey();
695  }
696 
697  }
698 
700  {
702  };
703 
704 
706  {
710  "BasisType is not a boundary interior form");
714  "BasisType is not a boundary interior form");
715 
716  return 4 + 2*(GetBasisNumModes(0)-2) + 2*(GetBasisNumModes(1)-2);
717  }
718 
720  {
724  "BasisType is not a boundary interior form");
728  "BasisType is not a boundary interior form");
729 
730  return 2*GetBasisNumModes(0) + 2*GetBasisNumModes(1);
731  }
732 
734  const std::vector<unsigned int> &nummodes,
735  int &modes_offset)
736  {
737  int nmodes = nummodes[modes_offset]*nummodes[modes_offset+1];
738  modes_offset += 2;
739 
740  return nmodes;
741  }
742 
744  {
745  bool returnval = false;
746 
749  {
752  {
753  returnval = true;
754  }
755  }
756 
757  return returnval;
758  }
759 
761  Array<OneD, NekDouble> &coords_1,
762  Array<OneD, NekDouble> &coords_2)
763  {
764  Array<OneD, const NekDouble> z0 = m_base[0]->GetZ();
765  Array<OneD, const NekDouble> z1 = m_base[1]->GetZ();
766  int nq0 = GetNumPoints(0);
767  int nq1 = GetNumPoints(1);
768  int i;
769 
770  for(i = 0; i < nq1; ++i)
771  {
772  Blas::Dcopy(nq0,z0.get(), 1,&coords_0[0] + i*nq0,1);
773  Vmath::Fill(nq0,z1[i],&coords_1[0] + i*nq0,1);
774  }
775  }
776 
777  //////////////
778  // Mappings //
779  //////////////
780 
782  {
783  int i;
784  int cnt=0;
785  int nummodes0, nummodes1;
786  int value1 = 0, value2 = 0;
787  if(outarray.num_elements()!=NumBndryCoeffs())
788  {
790  }
791 
792  nummodes0 = m_base[0]->GetNumModes();
793  nummodes1 = m_base[1]->GetNumModes();
794 
795  const LibUtilities::BasisType Btype0 = GetBasisType(0);
796  const LibUtilities::BasisType Btype1 = GetBasisType(1);
797 
798  switch(Btype1)
799  {
802  value1 = nummodes0;
803  break;
805  value1 = 2*nummodes0;
806  break;
807  default:
808  ASSERTL0(0,"Mapping array is not defined for this expansion");
809  break;
810  }
811 
812  for(i = 0; i < value1; i++)
813  {
814  outarray[i]=i;
815  }
816  cnt=value1;
817 
818  switch(Btype0)
819  {
822  value2 = value1+nummodes0-1;
823  break;
825  value2 = value1+1;
826  break;
827  default:
828  ASSERTL0(0,"Mapping array is not defined for this expansion");
829  break;
830  }
831 
832  for(i = 0; i < nummodes1-2; i++)
833  {
834  outarray[cnt++]=value1+i*nummodes0;
835  outarray[cnt++]=value2+i*nummodes0;
836  }
837 
838 
840  {
841  for(i = nummodes0*(nummodes1-1);i < GetNcoeffs(); i++)
842  {
843  outarray[cnt++] = i;
844  }
845  }
846  }
847 
849  {
850  int i,j;
851  int cnt=0;
852  int nummodes0, nummodes1;
853  int startvalue;
854  if(outarray.num_elements()!=GetNcoeffs()-NumBndryCoeffs())
855  {
857  }
858 
859  nummodes0 = m_base[0]->GetNumModes();
860  nummodes1 = m_base[1]->GetNumModes();
861 
862  const LibUtilities::BasisType Btype0 = GetBasisType(0);
863  const LibUtilities::BasisType Btype1 = GetBasisType(1);
864 
865  switch(Btype1)
866  {
868  startvalue = nummodes0;
869  break;
871  startvalue = 2*nummodes0;
872  break;
873  default:
874  ASSERTL0(0,"Mapping array is not defined for this expansion");
875  break;
876  }
877 
878  switch(Btype0)
879  {
881  startvalue++;
882  break;
884  startvalue+=2;
885  break;
886  default:
887  ASSERTL0(0,"Mapping array is not defined for this expansion");
888  break;
889  }
890 
891  for(i = 0; i < nummodes1-2; i++)
892  {
893  for(j = 0; j < nummodes0-2; j++)
894  {
895  outarray[cnt++]=startvalue+j;
896  }
897  startvalue+=nummodes0;
898  }
899  }
900 
901  int StdQuadExp::v_GetVertexMap(int localVertexId, bool useCoeffPacking)
902  {
903  int localDOF = 0;
904 
905  if(useCoeffPacking == true)
906  {
907  switch(localVertexId)
908  {
909  case 0:
910  {
911  localDOF = 0;
912  }
913  break;
914  case 1:
915  {
917  {
918  localDOF = m_base[0]->GetNumModes()-1;
919  }
920  else
921  {
922  localDOF = 1;
923  }
924  }
925  break;
926  case 2:
927  {
929  {
930  localDOF = m_base[0]->GetNumModes() * (m_base[1]->GetNumModes()-1);
931  }
932  else
933  {
934  localDOF = m_base[0]->GetNumModes();
935  }
936  }
937  break;
938  case 3:
939  {
941  {
942  localDOF = m_base[0]->GetNumModes()*m_base[1]->GetNumModes()-1;
943  }
944  else
945  {
946  localDOF = m_base[0]->GetNumModes()+1;
947  }
948  }
949  break;
950  default:
951  ASSERTL0(false,"eid must be between 0 and 3");
952  break;
953  }
954 
955  }
956  else
957  {
958  switch(localVertexId)
959  {
960  case 0:
961  {
962  localDOF = 0;
963  }
964  break;
965  case 1:
966  {
968  {
969  localDOF = m_base[0]->GetNumModes()-1;
970  }
971  else
972  {
973  localDOF = 1;
974  }
975  }
976  break;
977  case 2:
978  {
980  {
981  localDOF = m_base[0]->GetNumModes()*m_base[1]->GetNumModes()-1;
982  }
983  else
984  {
985  localDOF = m_base[0]->GetNumModes()+1;
986  }
987  }
988  break;
989  case 3:
990  {
992  {
993  localDOF = m_base[0]->GetNumModes() * (m_base[1]->GetNumModes()-1);
994  }
995  else
996  {
997  localDOF = m_base[0]->GetNumModes();
998  }
999  }
1000  break;
1001  default:
1002  ASSERTL0(false,"eid must be between 0 and 3");
1003  break;
1004  }
1005  }
1006  return localDOF;
1007  }
1008 
1010  const Orientation edgeOrient,
1011  Array<OneD, unsigned int> &maparray,
1012  Array<OneD, int> &signarray)
1013  {
1014  int i;
1015  const int nummodes0 = m_base[0]->GetNumModes();
1016  const int nummodes1 = m_base[1]->GetNumModes();
1017  const int nEdgeIntCoeffs = GetEdgeNcoeffs(eid)-2;
1018  const LibUtilities::BasisType bType = GetEdgeBasisType(eid);
1019 
1020  if(maparray.num_elements() != nEdgeIntCoeffs)
1021  {
1022  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1023  }
1024 
1025  if(signarray.num_elements() != nEdgeIntCoeffs)
1026  {
1027  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
1028  }
1029  else
1030  {
1031  fill( signarray.get() , signarray.get()+nEdgeIntCoeffs, 1 );
1032  }
1033 
1034  if(bType == LibUtilities::eModified_A)
1035  {
1036  switch(eid)
1037  {
1038  case 0:
1039  {
1040  for(i = 0; i < nEdgeIntCoeffs; i++)
1041  {
1042  maparray[i] = i+2;
1043  }
1044 
1045  if(edgeOrient==eBackwards)
1046  {
1047  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1048  {
1049  signarray[i] = -1;
1050  }
1051  }
1052  }
1053  break;
1054  case 1:
1055  {
1056  for(i = 0; i < nEdgeIntCoeffs; i++)
1057  {
1058  maparray[i] = (i+2)*nummodes0 + 1;
1059  }
1060 
1061  if(edgeOrient==eBackwards)
1062  {
1063  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1064  {
1065  signarray[i] = -1;
1066  }
1067  }
1068  }
1069  break;
1070  case 2:
1071  {
1072  for(i = 0; i < nEdgeIntCoeffs; i++)
1073  {
1074  maparray[i] = nummodes0+i+2;
1075  }
1076 
1077  if(edgeOrient==eForwards)
1078  {
1079  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1080  {
1081  signarray[i] = -1;
1082  }
1083  }
1084  }
1085  break;
1086  case 3:
1087  {
1088  for(i = 0; i < nEdgeIntCoeffs; i++)
1089  {
1090  maparray[i] = (i+2)*nummodes0;
1091  }
1092 
1093  if(edgeOrient==eForwards)
1094  {
1095  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1096  {
1097  signarray[i] = -1;
1098  }
1099  }
1100  }
1101  break;
1102  default:
1103  ASSERTL0(false,"eid must be between 0 and 3");
1104  break;
1105  }
1106  }
1107  else if(bType == LibUtilities::eGLL_Lagrange)
1108  {
1109  switch(eid)
1110  {
1111  case 0:
1112  {
1113  for(i = 0; i < nEdgeIntCoeffs; i++)
1114  {
1115  maparray[i] = i+1;
1116  }
1117  }
1118  break;
1119  case 1:
1120  {
1121  for(i = 0; i < nEdgeIntCoeffs; i++)
1122  {
1123  maparray[i] = (i+2)*nummodes0 - 1;
1124  }
1125  }
1126  break;
1127  case 2:
1128  {
1129  for(i = 0; i < nEdgeIntCoeffs; i++)
1130  {
1131  maparray[i] = nummodes0*nummodes1 - 2 - i;
1132  }
1133  }
1134  break;
1135  case 3:
1136  {
1137  for(i = 0; i < nEdgeIntCoeffs; i++)
1138  {
1139  maparray[i] = nummodes0*(nummodes1-2-i);
1140  }
1141  }
1142  break;
1143  default:
1144  ASSERTL0(false,"eid must be between 0 and 3");
1145  break;
1146  }
1147  if(edgeOrient == eBackwards)
1148  {
1149  reverse( maparray.get() , maparray.get()+nEdgeIntCoeffs );
1150  }
1151  }
1152  else
1153  {
1154  ASSERTL0(false,"Mapping not defined for this type of basis");
1155  }
1156 
1157  }
1158 
1160  const int eid,
1161  const Orientation edgeOrient,
1162  Array<OneD, unsigned int>& maparray,
1163  Array<OneD, int>& signarray,
1164  int P)
1165  {
1166  int i;
1167  int numModes=0;
1168  int order0 = m_base[0]->GetNumModes();
1169  int order1 = m_base[1]->GetNumModes();
1170 
1171  switch (eid)
1172  {
1173  case 0:
1174  case 2:
1175  numModes = order0;
1176  break;
1177  case 1:
1178  case 3:
1179  numModes = order1;
1180  }
1181 
1182  bool checkForZeroedModes = false;
1183  if (P == -1)
1184  {
1185  P = numModes;
1186  }
1187  else if(P != numModes)
1188  {
1189  checkForZeroedModes = true;
1190  }
1191  const LibUtilities::BasisType bType = GetEdgeBasisType(eid);
1192 
1193 
1194  if (maparray.num_elements() != P)
1195  {
1196  maparray = Array<OneD, unsigned int>(P);
1197  }
1198 
1199  if(signarray.num_elements() != P)
1200  {
1201  signarray = Array<OneD, int>(P, 1);
1202  }
1203  else
1204  {
1205  fill(signarray.get(), signarray.get()+P, 1);
1206  }
1207 
1208  if (bType == LibUtilities::eModified_A)
1209  {
1210  switch (eid)
1211  {
1212  case 0:
1213  {
1214  for (i = 0; i < P; i++)
1215  {
1216  maparray[i] = i;
1217  }
1218 
1219  if (edgeOrient == eBackwards)
1220  {
1221  swap(maparray[0], maparray[1]);
1222 
1223  for(i = 3; i < P; i+=2)
1224  {
1225  signarray[i] = -1;
1226  }
1227  }
1228  }
1229  break;
1230  case 1:
1231  {
1232  for (i = 0; i < P; i++)
1233  {
1234  maparray[i] = i*order0 + 1;
1235  }
1236 
1237  if (edgeOrient == eBackwards)
1238  {
1239  swap(maparray[0], maparray[1]);
1240 
1241  for(i = 3; i < P; i+=2)
1242  {
1243  signarray[i] = -1;
1244  }
1245  }
1246  }
1247  break;
1248  case 2:
1249  {
1250  for (i = 0; i < P; i++)
1251  {
1252  maparray[i] = order0+i;
1253  }
1254 
1255  if (edgeOrient == eForwards)
1256  {
1257  swap(maparray[0], maparray[1]);
1258 
1259  for (i = 3; i < P; i+=2)
1260  {
1261  signarray[i] = -1;
1262  }
1263  }
1264  }
1265  break;
1266  case 3:
1267  {
1268  for (i = 0; i < P; i++)
1269  {
1270  maparray[i] = i*order0;
1271  }
1272 
1273  if (edgeOrient == eForwards)
1274  {
1275  swap(maparray[0], maparray[1]);
1276 
1277  for (i = 3; i < P; i+=2)
1278  {
1279  signarray[i] = -1;
1280  }
1281  }
1282  }
1283  break;
1284  default:
1285  ASSERTL0(false, "eid must be between 0 and 3");
1286  break;
1287  }
1288  }
1289  else if(bType == LibUtilities::eGLL_Lagrange ||
1291  {
1292  switch (eid)
1293  {
1294  case 0:
1295  {
1296  for (i = 0; i < P; i++)
1297  {
1298  maparray[i] = i;
1299  }
1300  }
1301  break;
1302  case 1:
1303  {
1304  for (i = 0; i < P; i++)
1305  {
1306  maparray[i] = (i+1)*order0 - 1;
1307  }
1308  }
1309  break;
1310  case 2:
1311  {
1312  for (i = 0; i < P; i++)
1313  {
1314  maparray[i] = order0*order1 - 1 - i;
1315  }
1316  }
1317  break;
1318  case 3:
1319  {
1320  for (i = 0; i < P; i++)
1321  {
1322  maparray[i] = order0*(order1-1-i);
1323  }
1324  }
1325  break;
1326  default:
1327  ASSERTL0(false, "eid must be between 0 and 3");
1328  break;
1329  }
1330  if (edgeOrient == eBackwards)
1331  {
1332  reverse(maparray.get(), maparray.get()+P);
1333  }
1334  }
1335  else
1336  {
1337  ASSERTL0(false, "Mapping not defined for this type of basis");
1338  }
1339 
1340  if (checkForZeroedModes)
1341  {
1342  if (bType == LibUtilities::eModified_A)
1343  {
1344  // Zero signmap and set maparray to zero if
1345  // elemental modes are not as large as face modesl
1346  for (int j = numModes; j < P; j++)
1347  {
1348  signarray[j] = 0.0;
1349  maparray[j] = maparray[0];
1350  }
1351  }
1352  else
1353  {
1354  ASSERTL0(false, "Different trace space edge dimension "
1355  "and element edge dimension not possible "
1356  "for GLL-Lagrange bases");
1357  }
1358  }
1359  }
1360 
1361 
1362 
1363  ///////////////////////
1364  // Wrapper Functions //
1365  ///////////////////////
1366 
1368  {
1369  int i;
1370  int order0 = GetBasisNumModes(0);
1371  int order1 = GetBasisNumModes(1);
1372  MatrixType mtype = mkey.GetMatrixType();
1373 
1374  DNekMatSharedPtr Mat;
1375 
1376  switch(mtype)
1377  {
1379  {
1380  int nq0 = m_base[0]->GetNumPoints();
1381  int nq1 = m_base[1]->GetNumPoints();
1382  int nq;
1383 
1384  // take definition from key
1386  {
1387  nq = (int) mkey.GetConstFactor(eFactorConst);
1388  }
1389  else
1390  {
1391  nq = max(nq0,nq1);
1392  }
1393 
1394  int neq = LibUtilities::StdQuadData::
1395  getNumberOfCoefficients(nq, nq);
1396  Array<OneD, Array<OneD, NekDouble> > coords(neq);
1397  Array<OneD, NekDouble> coll (2);
1399  Array<OneD, NekDouble> tmp (nq0);
1400 
1402  AllocateSharedPtr(neq, nq0 * nq1);
1403  int cnt = 0;
1404 
1405  for(int i = 0; i < nq; ++i)
1406  {
1407  for(int j = 0; j < nq; ++j,++cnt)
1408  {
1409  coords[cnt] = Array<OneD, NekDouble>(2);
1410  coords[cnt][0] = -1.0 + 2*j/(NekDouble)(nq-1);
1411  coords[cnt][1] = -1.0 + 2*i/(NekDouble)(nq-1);
1412  }
1413  }
1414 
1415  for(int i = 0; i < neq; ++i)
1416  {
1417  LocCoordToLocCollapsed(coords[i],coll);
1418 
1419  I[0] = m_base[0]->GetI(coll);
1420  I[1] = m_base[1]->GetI(coll+1);
1421 
1422  // interpolate first coordinate direction
1423  for (int j = 0; j < nq1; ++j)
1424  {
1425  NekDouble fac = (I[1]->GetPtr())[j];
1426  Vmath::Smul(nq0,fac,I[0]->GetPtr(),1,tmp,1);
1427 
1428  Vmath::Vcopy(nq0, &tmp[0], 1,
1429  Mat->GetRawPtr()+j*nq0*neq+i,neq);
1430  }
1431 
1432  }
1433  break;
1434  }
1435  case eMass:
1436  {
1438  // For Fourier basis set the imaginary component of mean mode
1439  // to have a unit diagonal component in mass matrix
1441  {
1442  for(i = 0; i < order1; ++i)
1443  {
1444  (*Mat)(order0*i+1,i*order0+1) = 1.0;
1445  }
1446  }
1447 
1449  {
1450  for(i = 0; i < order0; ++i)
1451  {
1452  (*Mat)(order0+i ,order0+i) = 1.0;
1453  }
1454  }
1455  break;
1456  }
1457  case eFwdTrans:
1458  {
1460  StdMatrixKey iprodkey(eIProductWRTBase,DetShapeType(),*this);
1461  DNekMat &Iprod = *GetStdMatrix(iprodkey);
1462  StdMatrixKey imasskey(eInvMass,DetShapeType(),*this);
1463  DNekMat &Imass = *GetStdMatrix(imasskey);
1464 
1465  (*Mat) = Imass*Iprod;
1466  break;
1467  }
1468  case eGaussDG:
1469  {
1470  ConstFactorMap factors = mkey.GetConstFactors();
1471 
1472  int edge = (int)factors[StdRegions::eFactorGaussEdge];
1473  int dir = (edge + 1) % 2;
1474  int nCoeffs = m_base[dir]->GetNumModes();
1475 
1476  const LibUtilities::PointsKey BS_p(
1478  const LibUtilities::BasisKey BS_k(
1479  LibUtilities::eGauss_Lagrange, nCoeffs, BS_p);
1480 
1481  Array<OneD, NekDouble> coords(1, 0.0);
1482  coords[0] = (edge == 0 || edge == 3) ? -1.0 : 1.0;
1483 
1486  DNekMatSharedPtr m_Ix = basis->GetI(coords);
1487 
1489  1.0, nCoeffs);
1490  Vmath::Vcopy(nCoeffs, m_Ix->GetPtr(), 1, Mat->GetPtr(), 1);
1491  break;
1492  }
1493  default:
1494  {
1496  break;
1497  }
1498  }
1499 
1500  return Mat;
1501  }
1502 
1504  {
1505  return GenMatrix(mkey);
1506  }
1507 
1508  ///////////////////////////////////
1509  // Operator evaluation functions //
1510  ///////////////////////////////////
1511 
1513  const Array<OneD, const NekDouble> &inarray,
1514  Array<OneD,NekDouble> &outarray,
1515  const StdMatrixKey &mkey)
1516  {
1517 
1519 
1520  if(inarray.get() == outarray.get())
1521  {
1523  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
1524 
1525  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, 1.0, mat->GetPtr().get(),
1526  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
1527  }
1528  else
1529  {
1530  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, 1.0, mat->GetPtr().get(),
1531  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
1532  }
1533  }
1534 
1535 
1537  const StdMatrixKey &mkey)
1538  {
1539  // Generate an orthonogal expansion
1540  int qa = m_base[0]->GetNumPoints();
1541  int qb = m_base[1]->GetNumPoints();
1542  int nmodes_a = m_base[0]->GetNumModes();
1543  int nmodes_b = m_base[1]->GetNumModes();
1544  int nmodes = min(nmodes_a,nmodes_b);
1545  // Declare orthogonal basis.
1548 
1551  StdQuadExp OrthoExp(Ba,Bb);
1552 
1553  //SVV parameters loaded from the .xml case file
1554  NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDiffCoeff);
1555  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
1556 
1557  if(mkey.HasVarCoeff(eVarCoeffLaplacian)) // Rodrigo's svv mapping
1558  {
1559  Array<OneD, NekDouble> sqrt_varcoeff(qa*qb);
1560  Array<OneD, NekDouble> tmp(qa*qb);
1561 
1562  Vmath::Vsqrt(qa * qb,
1564  sqrt_varcoeff, 1);
1565 
1566  //Vmath::Fill(qa*qb,Vmath::Vmax(qa*qb,sqrt_varcoeff,1),
1567  //sqrt_varcoeff,1);
1568 
1569  // multiply by sqrt(Variable Coefficient) containing h v /p
1570  Vmath::Vmul(qa*qb,sqrt_varcoeff,1,array,1,tmp,1);
1571 
1572  // project onto modal space.
1573  OrthoExp.FwdTrans(tmp,orthocoeffs);
1574 
1575  for(int j = 0; j < nmodes_a; ++j)
1576  {
1577  for(int k = 0; k < nmodes_b; ++k)
1578  {
1579  // linear space but makes high modes very negative
1580  orthocoeffs[j*nmodes_b+k] *=
1581  (1.0+SvvDiffCoeff*
1582  pow(j/(nmodes_a-1)+k/(nmodes_b-1),0.5*nmodes));
1583  // bilinear blend
1584  //orthocoeffs[j*nmodes_b+k] *=
1585  //(1.0 + SvvDiffCoeff
1586  // *pow(j/(NekDouble)(nmodes_a-1.0),0.5*nmodes)
1587  // *pow(k/(NekDouble)(nmodes_b-1.0),0.5*nmodes));
1588  }
1589  }
1590 
1591  // backward transform to physical space
1592  OrthoExp.BwdTrans(orthocoeffs,tmp);
1593 
1594  // multiply by sqrt(Variable Coefficient) containing h v /p
1595  // - split to keep symmetry
1596  Vmath::Vmul(qa*qb,sqrt_varcoeff,1,tmp,1,array,1);
1597  }
1598  else
1599  {
1600  //for the "old" implementation
1601  int cutoff = (int) (mkey.GetConstFactor(eFactorSVVCutoffRatio)*
1602  min(nmodes_a,nmodes_b));
1603 
1604  // project onto modal space.
1605  OrthoExp.FwdTrans(array,orthocoeffs);
1606 
1607  //counters for scanning through orthocoeffs array
1608  int j, k, cnt = 0;
1609 
1610  //------"New" Version August 22nd '13--------------------
1611  for(j = 0; j < nmodes_a; ++j)
1612  {
1613  for(k = 0; k < nmodes_b; ++k)
1614  {
1615  if(j + k >= cutoff)//to filter out only the "high-modes"
1616  {
1617  orthocoeffs[j*nmodes_b+k] *=
1618  (SvvDiffCoeff*exp(-(j+k-nmodes)*(j+k-nmodes)/
1619  ((NekDouble)((j+k-cutoff+1)*
1620  (j+k-cutoff+1)))));
1621  }
1622  else
1623  {
1624  orthocoeffs[j*nmodes_b+k] *= 0.0;
1625  }
1626  cnt++;
1627  }
1628  }
1629 
1630  // backward transform to physical space
1631  OrthoExp.BwdTrans(orthocoeffs,array);
1632  }
1633  }
1634 
1636  int numMin,
1637  const Array<OneD, const NekDouble> &inarray,
1638  Array<OneD, NekDouble> &outarray)
1639  {
1640  int n_coeffs = inarray.num_elements();
1641 
1642 
1643  Array<OneD, NekDouble> coeff(n_coeffs);
1644  Array<OneD, NekDouble> coeff_tmp(n_coeffs,0.0);
1647 
1648  int nmodes0 = m_base[0]->GetNumModes();
1649  int nmodes1 = m_base[1]->GetNumModes();
1650  int numMax = nmodes0;
1651 
1652  Vmath::Vcopy(n_coeffs,inarray,1,coeff_tmp,1);
1653 
1654  const LibUtilities::PointsKey Pkey0(
1656  const LibUtilities::PointsKey Pkey1(
1658 
1659  LibUtilities::BasisKey b0(m_base[0]->GetBasisType(),nmodes0,Pkey0);
1660  LibUtilities::BasisKey b1(m_base[1]->GetBasisType(),nmodes1,Pkey1);
1661 
1662  LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A,nmodes0,Pkey0);
1663  LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_A,nmodes1,Pkey1);
1664 
1666  b0, b1, coeff_tmp, bortho0, bortho1, coeff);
1667 
1668  Vmath::Zero(n_coeffs,coeff_tmp,1);
1669 
1670  int cnt = 0;
1671  for (int i = 0; i < numMin+1; ++i)
1672  {
1673  Vmath::Vcopy(numMin,
1674  tmp = coeff+cnt,1,
1675  tmp2 = coeff_tmp+cnt,1);
1676 
1677  cnt = i*numMax;
1678  }
1679 
1681  bortho0, bortho1, coeff_tmp, b0, b1, outarray);
1682  }
1683 
1684 
1686  const Array<OneD, const NekDouble> &inarray,
1687  Array<OneD,NekDouble> &outarray,
1688  const StdMatrixKey &mkey)
1689  {
1690  StdExpansion::MassMatrixOp_MatFree(inarray,outarray,mkey);
1691  }
1692 
1694  const Array<OneD, const NekDouble> &inarray,
1695  Array<OneD,NekDouble> &outarray,
1696  const StdMatrixKey &mkey)
1697  {
1698  StdQuadExp::v_LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
1699  }
1700 
1701  void StdQuadExp::v_LaplacianMatrixOp(const int k1, const int k2,
1702  const Array<OneD, const NekDouble> &inarray,
1703  Array<OneD,NekDouble> &outarray,
1704  const StdMatrixKey &mkey)
1705  {
1706  StdExpansion::LaplacianMatrixOp_MatFree(k1,k2,inarray,outarray,mkey);
1707  }
1708 
1710  const Array<OneD, const NekDouble> &inarray,
1711  Array<OneD,NekDouble> &outarray,
1712  const StdMatrixKey &mkey)
1713  {
1714  StdExpansion::WeakDerivMatrixOp_MatFree(i,inarray,outarray,mkey);
1715  }
1716 
1718  Array<OneD,NekDouble> &outarray,
1719  const StdMatrixKey &mkey)
1720  {
1721  StdQuadExp::v_HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
1722  }
1723 
1724  //up to here
1726  const Array<OneD, const NekDouble> &inarray,
1727  Array<OneD, NekDouble> &outarray)
1728  {
1729  int i;
1730  int nquad0 = m_base[0]->GetNumPoints();
1731  int nquad1 = m_base[1]->GetNumPoints();
1732 
1733  const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
1734  const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
1735 
1736  // multiply by integration constants
1737  for(i = 0; i < nquad1; ++i)
1738  {
1739  Vmath::Vmul(nquad0, inarray.get()+i*nquad0,1,
1740  w0.get(),1,outarray.get()+i*nquad0,1);
1741  }
1742 
1743  for(i = 0; i < nquad0; ++i)
1744  {
1745  Vmath::Vmul(nquad1,outarray.get()+i,nquad0,w1.get(),1,
1746  outarray.get()+i,nquad0);
1747  }
1748  }
1749 
1751  Array<OneD, int> &conn,
1752  bool standard)
1753  {
1754  int np1 = m_base[0]->GetNumPoints();
1755  int np2 = m_base[1]->GetNumPoints();
1756  int np = max(np1,np2);
1757 
1758  conn = Array<OneD, int>(6*(np-1)*(np-1));
1759 
1760  int row = 0;
1761  int rowp1 = 0;
1762  int cnt = 0;
1763  for(int i = 0; i < np-1; ++i)
1764  {
1765  rowp1 += np;
1766  for(int j = 0; j < np-1; ++j)
1767  {
1768  conn[cnt++] = row +j;
1769  conn[cnt++] = row +j+1;
1770  conn[cnt++] = rowp1 +j;
1771 
1772  conn[cnt++] = rowp1 +j+1;
1773  conn[cnt++] = rowp1 +j;
1774  conn[cnt++] = row +j+1;
1775  }
1776  row += np;
1777  }
1778  }
1779 
1780  } //end namespace
1781 }//end namespace
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:470
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:122
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
#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
const ConstFactorMap & GetConstFactors() const
Definition: StdMatrixKey.h:142
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdQuadExp.cpp:443
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:82
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
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: StdQuadExp.cpp:97
void IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
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 DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
this function calculates the inner product of a given function f with the different modes of the expa...
Definition: StdExpansion.h:635
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
const Array< OneD, const NekDouble > & GetVarCoeff(const StdRegions::VarCoeffType &coeff) const
Definition: StdMatrixKey.h:152
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
Definition: StdExpansion.h:229
Principle Modified Functions .
Definition: BasisType.h:49
virtual bool v_IsBoundaryInteriorExpansion()
Definition: StdQuadExp.cpp:743
Lagrange Polynomials using the Gauss points .
Definition: BasisType.h:54
STL namespace.
void LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
Convert local cartesian coordinate xi into local collapsed coordinates eta.
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:252
NekDouble Integral(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &w0, const Array< OneD, const NekDouble > &w1)
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdQuadExp.cpp:848
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
virtual void v_LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Fourier Expansion .
Definition: BasisType.h:52
LibUtilities::NekManager< StdMatrixKey, DNekBlkMat, StdMatrixKey::opLess > m_stdStaticCondMatrixManager
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_StdPhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
Definition: StdQuadExp.cpp:129
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...
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:49
virtual void v_FwdTrans_BndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdQuadExp.cpp:256
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: StdQuadExp.cpp:233
void v_GeneralMatrixOp_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:131
BasisManagerT & BasisManager(void)
virtual int v_NumDGBndryCoeffs() const
Definition: StdQuadExp.cpp:719
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:141
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdQuadExp.cpp:781
LibUtilities::BasisType GetEdgeBasisType(const int i) const
This function returns the type of expansion basis on the i-th edge.
Definition: StdExpansion.h:413
void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
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: StdQuadExp.cpp:509
The base class for all shapes.
Definition: StdExpansion.h:69
virtual const LibUtilities::BasisKey v_DetEdgeBasisKey(const int i) const
Definition: StdQuadExp.cpp:684
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_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
virtual void v_IProductWRTBase_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdQuadExp.cpp:424
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdQuadExp.cpp:166
virtual int v_DetCartesianDirOfEdge(const int edge)
Definition: StdQuadExp.cpp:670
void WeakDerivMatrixOp_MatFree(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
void GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:819
virtual void v_GetEdgeInteriorMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray)
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
Principle Orthogonal Functions .
Definition: BasisType.h:46
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdQuadExp.cpp:436
void v_GetEdgeToElementMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, int P=-1)
virtual void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
Definition: StdQuadExp.cpp:569
Defines a specification for a set of points.
Definition: Points.h:58
virtual void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
double NekDouble
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_0, Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2)
Definition: StdQuadExp.cpp:760
virtual int v_GetEdgeNcoeffs(const int i) const
Definition: StdQuadExp.cpp:628
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 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: StdQuadExp.cpp:188
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_IProductWRTDerivBase_MatOp(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdQuadExp.cpp:474
virtual int v_GetNedges() const
Definition: StdQuadExp.cpp:623
virtual int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset)
Definition: StdQuadExp.cpp:733
void GetEdgeToElementMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, int P=-1)
Definition: StdExpansion.h:851
virtual LibUtilities::BasisType v_GetEdgeBasisType(const int i) const
Definition: StdQuadExp.cpp:656
virtual int v_NumBndryCoeffs() const
Definition: StdQuadExp.cpp:705
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:132
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:250
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*base1 and put into outarray...
Definition: StdQuadExp.cpp:381
void MassMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
Definition: StdQuadExp.cpp:901
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 void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Definition: StdQuadExp.cpp:395
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
NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
Integrates the specified function over the domain.
Definition: StdQuadExp.cpp:78
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdQuadExp.cpp:152
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &array)
Fill outarray with mode mode of expansion.
Definition: StdQuadExp.cpp:582
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)
virtual void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
boost::shared_ptr< Basis > BasisSharedPtr
Lagrange for SEM basis .
Definition: BasisType.h:53
virtual int v_GetEdgeNumPoints(const int i) const
Definition: StdQuadExp.cpp:642
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
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 InterpCoeff2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
Definition: InterpCoeff.cpp:73
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
Describes the specification for a Basis.
Definition: Basis.h:50
virtual LibUtilities::ShapeType v_DetShapeType() const
Definition: StdQuadExp.cpp:699
virtual int v_GetNverts() const
Definition: StdQuadExp.cpp:618
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:52
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_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true)
virtual void v_HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
virtual void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space...