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