Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 ///////////////////////////////////////////////////////////////////////////////
35 #include <LibUtilities/Foundations/InterpCoeff.h>
36 #include <StdRegions/StdQuadExp.h>
37 #include <StdRegions/StdSegExp.h>
38 #include <LibUtilities/Foundations/ManagerAccess.h>
39 
40 namespace Nektar
41 {
42  namespace StdRegions
43  {
44 
45 
46  StdQuadExp::StdQuadExp()
47  {
48  }
49 
50  /** \brief Constructor using BasisKey class for quadrature
51  * points and order definition
52  */
53  StdQuadExp::StdQuadExp(const LibUtilities::BasisKey &Ba,
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 */
61  StdQuadExp::StdQuadExp(const StdQuadExp &T):
62  StdExpansion(T),
63  StdExpansion2D(T)
64  {
65  }
66 
67  /** \brief Destructor */
68  StdQuadExp::~StdQuadExp()
69  {
70  }
71 
72  /////////////////////////
73  // Integration Methods //
74  /////////////////////////
75 
76  NekDouble StdQuadExp::v_Integral(
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 
95  void StdQuadExp::v_PhysDeriv(const Array<OneD, const NekDouble>& inarray,
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 
127  void StdQuadExp::v_StdPhysDeriv(const Array<OneD, const NekDouble>& inarray,
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 
150  void StdQuadExp::v_BwdTrans(const Array<OneD, const NekDouble>& inarray,
151  Array<OneD, NekDouble> &outarray)
152  {
153  if(m_base[0]->Collocation() && m_base[1]->Collocation())
154  {
155  Vmath::Vcopy(m_base[0]->GetNumPoints() * m_base[1]->GetNumPoints(),
156  inarray, 1, outarray, 1);
157  }
158  else
159  {
160  StdQuadExp::v_BwdTrans_SumFac(inarray,outarray);
161  }
162  }
163 
164  void StdQuadExp::v_BwdTrans_SumFac(
165  const Array<OneD, const NekDouble>& inarray,
166  Array<OneD, NekDouble> &outarray)
167  {
168  Array<OneD, NekDouble> wsp(m_base[0]->GetNumPoints()*
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;
186  void StdQuadExp::v_BwdTrans_SumFacKernel(
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,
191  Array<OneD, NekDouble> &wsp,
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 
231  void StdQuadExp::v_FwdTrans(const Array<OneD, const NekDouble>& inarray,
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);
248  NekVector<NekDouble> out(m_ncoeffs,outarray,eWrapper);
249 
250  out = (*matsys)*in;
251  }
252  }
253 
254  void StdQuadExp::v_FwdTrans_BndConstrained(
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 
292  StdSegExpSharedPtr segexp[2] = {MemoryManager<StdRegions::StdSegExp>::AllocateSharedPtr(m_base[0]->GetBasisKey()),
293  MemoryManager<StdRegions::StdSegExp>::AllocateSharedPtr(m_base[1]->GetBasisKey())};
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 
311  Array<OneD, NekDouble> tmp0(m_ncoeffs);
312  Array<OneD, NekDouble> tmp1(m_ncoeffs);
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  */
379  void StdQuadExp::v_IProductWRTBase(
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 
393 
394  void StdQuadExp::v_IProductWRTBase_SumFac(
395  const Array<OneD, const NekDouble>& inarray,
396  Array<OneD, NekDouble> &outarray)
397  {
398  int nquad0 = m_base[0]->GetNumPoints();
399  int nquad1 = m_base[1]->GetNumPoints();
400  int order0 = m_base[0]->GetNumModes();
401 
402  Array<OneD,NekDouble> tmp(nquad0*nquad1+nquad1*order0);
403  Array<OneD,NekDouble> wsp(tmp+nquad0*nquad1);
404 
405  // multiply by integration constants
406  MultiplyByQuadratureMetric(inarray,tmp);
407 
408  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
409  m_base[1]->GetBdata(),
410  tmp,outarray,wsp,true,true);
411  }
412 
413  void StdQuadExp::v_IProductWRTBase_MatOp(
414  const Array<OneD, const NekDouble>& inarray,
415  Array<OneD, NekDouble> &outarray)
416  {
417  int nq = GetTotPoints();
418  StdMatrixKey iprodmatkey(eIProductWRTBase,DetShapeType(),*this);
419  DNekMatSharedPtr iprodmat = GetStdMatrix(iprodmatkey);
420 
421  Blas::Dgemv('N',m_ncoeffs,nq,1.0,iprodmat->GetPtr().get(),
422  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
423  }
424 
425  void StdQuadExp::v_IProductWRTDerivBase(const int dir,
426  const Array<OneD, const NekDouble>& inarray,
427  Array<OneD, NekDouble> & outarray)
428  {
429  StdQuadExp::IProductWRTDerivBase_SumFac(dir,inarray,outarray);
430  }
431 
432  void StdQuadExp::v_IProductWRTDerivBase_SumFac(const int dir,
433  const Array<OneD, const NekDouble>& inarray,
434  Array<OneD, NekDouble> &outarray)
435  {
436  ASSERTL0((dir==0)||(dir==1),"input dir is out of range");
437 
438  int nquad0 = m_base[0]->GetNumPoints();
439  int nquad1 = m_base[1]->GetNumPoints();
440  int nqtot = nquad0*nquad1;
441  int order0 = m_base[0]->GetNumModes();
442 
443  Array<OneD,NekDouble> tmp(nqtot+nquad1*order0);
444  Array<OneD,NekDouble> wsp(tmp+nqtot);
445 
446  // multiply by integration constants
447  MultiplyByQuadratureMetric(inarray,tmp);
448 
449  if(dir) // dir == 1
450  {
451  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
452  m_base[1]->GetDbdata(),
453  tmp,outarray,wsp,true,false);
454  }
455  else // dir == 0
456  {
457  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
458  m_base[1]->GetBdata(),
459  tmp,outarray,wsp,false,true);
460  }
461  }
462 
463  void StdQuadExp::v_IProductWRTDerivBase_MatOp(const int dir,
464  const Array<OneD, const NekDouble>& inarray,
465  Array<OneD, NekDouble> &outarray)
466  {
467  ASSERTL0((dir==0)||(dir==1),"input dir is out of range");
468 
469  int nq = GetTotPoints();
470  MatrixType mtype;
471 
472  if(dir) // dir == 1
473  {
474  mtype = eIProductWRTDerivBase1;
475  }
476  else // dir == 0
477  {
478  mtype = eIProductWRTDerivBase0;
479  }
480 
481  StdMatrixKey iprodmatkey(mtype,DetShapeType(),*this);
482  DNekMatSharedPtr iprodmat = GetStdMatrix(iprodmatkey);
483 
484  Blas::Dgemv('N',m_ncoeffs,nq,1.0,iprodmat->GetPtr().get(),
485  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
486  }
487 
488  // the arguments doCheckCollDir0 and doCheckCollDir1 allow you to specify whether
489  // to check if the basis has collocation properties (i.e. for the classical spectral
490  // element basis, In this case the 1D 'B' matrix is equal to the identity matrix
491  // which can be exploited to speed up the calculations).
492  // However, as this routine also allows to pass the matrix 'DB' (derivative of the basis),
493  // the collocation property cannot always be used. Therefor follow this rule:
494  // if base0 == m_base[0]->GetBdata() --> set doCheckCollDir0 == true;
495  // base1 == m_base[1]->GetBdata() --> set doCheckCollDir1 == true;
496  // base0 == m_base[0]->GetDbdata() --> set doCheckCollDir0 == false;
497  // base1 == m_base[1]->GetDbdata() --> set doCheckCollDir1 == false;
498  void StdQuadExp::v_IProductWRTBase_SumFacKernel(
499  const Array<OneD, const NekDouble>& base0,
500  const Array<OneD, const NekDouble>& base1,
501  const Array<OneD, const NekDouble>& inarray,
502  Array<OneD, NekDouble> &outarray,
503  Array<OneD, NekDouble> &wsp,
504  bool doCheckCollDir0,
505  bool doCheckCollDir1)
506  {
507  int nquad0 = m_base[0]->GetNumPoints();
508  int nquad1 = m_base[1]->GetNumPoints();
509  int nmodes0 = m_base[0]->GetNumModes();
510  int nmodes1 = m_base[1]->GetNumModes();
511 
512  bool colldir0 = doCheckCollDir0?(m_base[0]->Collocation()):false;
513  bool colldir1 = doCheckCollDir1?(m_base[1]->Collocation()):false;
514 
515  if(colldir0 && colldir1)
516  {
517  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,outarray.get(),1);
518  }
519  else if(colldir0)
520  {
521  Blas::Dgemm('N', 'N',nmodes0,nmodes1, nquad1,1.0, inarray.get(),
522  nmodes0, base1.get(), nquad1, 0.0,outarray.get(),nmodes0);
523  }
524  else if(colldir1)
525  {
526  Blas::Dgemm('T','N',nmodes0,nquad1,nquad0,1.0,base0.get(),
527  nquad0,inarray.get(),nquad0,0.0,outarray.get(),nmodes0);
528  }
529  else
530  {
531  ASSERTL1(wsp.num_elements()>=nquad1*nmodes0,"Workspace size is not sufficient");
532 
533 #if 1
534  Blas::Dgemm('T','N',nmodes0,nquad1,nquad0,1.0,base0.get(),
535  nquad0,inarray.get(),nquad0,0.0,wsp.get(),nmodes0);
536 
537 #else
538  for(int i = 0; i < nmodes0; ++i)
539  {
540  for(int j = 0; j < nquad1; ++j)
541  {
542  wsp[j*nmodes0+i] = Blas::Ddot(nquad0,
543  base0.get()+i*nquad0,1,
544  inarray.get()+j*nquad0,1);
545  }
546  }
547 #endif
548  Blas::Dgemm('N', 'N',nmodes0,nmodes1, nquad1,1.0, wsp.get(),
549  nmodes0, base1.get(), nquad1, 0.0,outarray.get(),nmodes0);
550  }
551  }
552 
553  //////////////////////////
554  // Evaluation functions //
555  //////////////////////////
556 
557 
558  void StdQuadExp::v_LocCoordToLocCollapsed(const Array<OneD, const NekDouble>& xi,
559  Array<OneD, NekDouble>& eta)
560  {
561  eta[0] = xi[0];
562  eta[1] = xi[1];
563  }
564 
565  /** \brief Fill outarray with mode \a mode of expansion
566  *
567  * Note for quadrilateral expansions _base[0] (i.e. p) modes run
568  * fastest
569  */
570 
571  void StdQuadExp::v_FillMode(const int mode,
572  Array<OneD, NekDouble> &outarray)
573  {
574  int i;
575  int nquad0 = m_base[0]->GetNumPoints();
576  int nquad1 = m_base[1]->GetNumPoints();
577  Array<OneD, const NekDouble> base0 = m_base[0]->GetBdata();
578  Array<OneD, const NekDouble> base1 = m_base[1]->GetBdata();
579  int btmp0 = m_base[0]->GetNumModes();
580  int mode0 = mode%btmp0;
581  int mode1 = mode/btmp0;
582 
583 
584  ASSERTL2(mode1 == (int)floor((1.0*mode)/btmp0),
585  "Integer Truncation not Equiv to Floor");
586 
587  ASSERTL2(m_ncoeffs <= mode,
588  "calling argument mode is larger than total expansion order");
589 
590  for(i = 0; i < nquad1; ++i)
591  {
592  Vmath::Vcopy(nquad0,(NekDouble *)(base0.get() + mode0*nquad0),
593  1, &outarray[0]+i*nquad0,1);
594  }
595 
596  for(i = 0; i < nquad0; ++i)
597  {
598  Vmath::Vmul(nquad1,(NekDouble *)(base1.get() + mode1*nquad1),1,
599  &outarray[0]+i,nquad0,&outarray[0]+i,nquad0);
600  }
601  }
602 
603  //////////////////////
604  // Helper functions //
605  //////////////////////
606 
607  int StdQuadExp::v_GetNverts() const
608  {
609  return 4;
610  }
611 
612  int StdQuadExp::v_GetNedges() const
613  {
614  return 4;
615  }
616 
617  int StdQuadExp::v_GetEdgeNcoeffs(const int i) const
618  {
619  ASSERTL2((i >= 0)&&(i <= 3),"edge id is out of range");
620 
621  if((i == 0)||(i == 2))
622  {
623  return GetBasisNumModes(0);
624  }
625  else
626  {
627  return GetBasisNumModes(1);
628  }
629  }
630 
631  int StdQuadExp::v_GetEdgeNumPoints(const int i) const
632  {
633  ASSERTL2((i >= 0)&&(i <= 3),"edge id is out of range");
634 
635  if((i == 0)||(i == 2))
636  {
637  return GetNumPoints(0);
638  }
639  else
640  {
641  return GetNumPoints(1);
642  }
643  }
644 
645  LibUtilities::BasisType StdQuadExp::v_GetEdgeBasisType(const int i) const
646  {
647  ASSERTL2((i >= 0)&&(i <= 3),"edge id is out of range");
648 
649  if((i == 0)||(i == 2))
650  {
651  return GetBasisType(0);
652  }
653  else
654  {
655  return GetBasisType(1);
656  }
657  }
658 
659  int StdQuadExp::v_DetCartesianDirOfEdge(const int edge)
660  {
661  ASSERTL2((edge >= 0)&&(edge <= 3),"edge id is out of range");
662 
663  if((edge == 0)||(edge == 2))
664  {
665  return 0;
666  }
667  else
668  {
669  return 1;
670  }
671  }
672 
673  const LibUtilities::BasisKey StdQuadExp::v_DetEdgeBasisKey(const int i) const
674  {
675  ASSERTL2((i >= 0)&&(i <= 3),"edge id is out of range");
676 
677  if((i == 0)||(i == 2))
678  {
679  return GetBasis(0)->GetBasisKey();
680  }
681  else
682  {
683  return GetBasis(1)->GetBasisKey();
684  }
685 
686  }
687 
688  LibUtilities::ShapeType StdQuadExp::v_DetShapeType() const
689  {
690  return LibUtilities::eQuadrilateral;
691  };
692 
693 
694  int StdQuadExp::v_NumBndryCoeffs() const
695  {
696  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
697  GetBasisType(0) == LibUtilities::eGLL_Lagrange||
698  GetBasisType(0) == LibUtilities::eGauss_Lagrange,
699  "BasisType is not a boundary interior form");
700  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
701  GetBasisType(1) == LibUtilities::eGLL_Lagrange||
702  GetBasisType(0) == LibUtilities::eGauss_Lagrange,
703  "BasisType is not a boundary interior form");
704 
705  return 4 + 2*(GetBasisNumModes(0)-2) + 2*(GetBasisNumModes(1)-2);
706  }
707 
708  int StdQuadExp::v_NumDGBndryCoeffs() const
709  {
710  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
711  GetBasisType(0) == LibUtilities::eGLL_Lagrange ||
712  GetBasisType(0) == LibUtilities::eGauss_Lagrange,
713  "BasisType is not a boundary interior form");
714  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
715  GetBasisType(1) == LibUtilities::eGLL_Lagrange ||
716  GetBasisType(0) == LibUtilities::eGauss_Lagrange,
717  "BasisType is not a boundary interior form");
718 
719  return 2*GetBasisNumModes(0) + 2*GetBasisNumModes(1);
720  }
721 
722  int StdQuadExp::v_CalcNumberOfCoefficients(
723  const std::vector<unsigned int> &nummodes,
724  int &modes_offset)
725  {
726  int nmodes = nummodes[modes_offset]*nummodes[modes_offset+1];
727  modes_offset += 2;
728 
729  return nmodes;
730  }
731 
732  bool StdQuadExp::v_IsBoundaryInteriorExpansion()
733  {
734  bool returnval = false;
735 
736  if((m_base[0]->GetBasisType() == LibUtilities::eModified_A)
737  ||(m_base[0]->GetBasisType() == LibUtilities::eGLL_Lagrange))
738  {
739  if((m_base[1]->GetBasisType() == LibUtilities::eModified_A)
740  ||(m_base[1]->GetBasisType() == LibUtilities::eGLL_Lagrange))
741  {
742  returnval = true;
743  }
744  }
745 
746  return returnval;
747  }
748 
749  void StdQuadExp::v_GetCoords(Array<OneD, NekDouble> &coords_0,
750  Array<OneD, NekDouble> &coords_1,
751  Array<OneD, NekDouble> &coords_2)
752  {
753  Array<OneD, const NekDouble> z0 = m_base[0]->GetZ();
754  Array<OneD, const NekDouble> z1 = m_base[1]->GetZ();
755  int nq0 = GetNumPoints(0);
756  int nq1 = GetNumPoints(1);
757  int i;
758 
759  for(i = 0; i < nq1; ++i)
760  {
761  Blas::Dcopy(nq0,z0.get(), 1,&coords_0[0] + i*nq0,1);
762  Vmath::Fill(nq0,z1[i],&coords_1[0] + i*nq0,1);
763  }
764  }
765 
766  //////////////
767  // Mappings //
768  //////////////
769 
770  void StdQuadExp::v_GetBoundaryMap(Array<OneD, unsigned int>& outarray)
771  {
772  int i;
773  int cnt=0;
774  int nummodes0, nummodes1;
775  int value1 = 0, value2 = 0;
776  if(outarray.num_elements()!=NumBndryCoeffs())
777  {
778  outarray = Array<OneD, unsigned int>(NumBndryCoeffs());
779  }
780 
781  nummodes0 = m_base[0]->GetNumModes();
782  nummodes1 = m_base[1]->GetNumModes();
783 
784  const LibUtilities::BasisType Btype0 = GetBasisType(0);
785  const LibUtilities::BasisType Btype1 = GetBasisType(1);
786 
787  switch(Btype1)
788  {
789  case LibUtilities::eGLL_Lagrange:
790  case LibUtilities::eGauss_Lagrange:
791  value1 = nummodes0;
792  break;
793  case LibUtilities::eModified_A:
794  value1 = 2*nummodes0;
795  break;
796  default:
797  ASSERTL0(0,"Mapping array is not defined for this expansion");
798  break;
799  }
800 
801  for(i = 0; i < value1; i++)
802  {
803  outarray[i]=i;
804  }
805  cnt=value1;
806 
807  switch(Btype0)
808  {
809  case LibUtilities::eGLL_Lagrange:
810  case LibUtilities::eGauss_Lagrange:
811  value2 = value1+nummodes0-1;
812  break;
813  case LibUtilities::eModified_A:
814  value2 = value1+1;
815  break;
816  default:
817  ASSERTL0(0,"Mapping array is not defined for this expansion");
818  break;
819  }
820 
821  for(i = 0; i < nummodes1-2; i++)
822  {
823  outarray[cnt++]=value1+i*nummodes0;
824  outarray[cnt++]=value2+i*nummodes0;
825  }
826 
827 
828  if(Btype1 == LibUtilities::eGLL_Lagrange || Btype1 == LibUtilities::eGauss_Lagrange )
829  {
830  for(i = nummodes0*(nummodes1-1);i < GetNcoeffs(); i++)
831  {
832  outarray[cnt++] = i;
833  }
834  }
835  }
836 
837  void StdQuadExp::v_GetInteriorMap(Array<OneD, unsigned int>& outarray)
838  {
839  int i,j;
840  int cnt=0;
841  int nummodes0, nummodes1;
842  int startvalue;
843  if(outarray.num_elements()!=GetNcoeffs()-NumBndryCoeffs())
844  {
845  outarray = Array<OneD, unsigned int>(GetNcoeffs()-NumBndryCoeffs());
846  }
847 
848  nummodes0 = m_base[0]->GetNumModes();
849  nummodes1 = m_base[1]->GetNumModes();
850 
851  const LibUtilities::BasisType Btype0 = GetBasisType(0);
852  const LibUtilities::BasisType Btype1 = GetBasisType(1);
853 
854  switch(Btype1)
855  {
856  case LibUtilities::eGLL_Lagrange:
857  startvalue = nummodes0;
858  break;
859  case LibUtilities::eModified_A:
860  startvalue = 2*nummodes0;
861  break;
862  default:
863  ASSERTL0(0,"Mapping array is not defined for this expansion");
864  break;
865  }
866 
867  switch(Btype0)
868  {
869  case LibUtilities::eGLL_Lagrange:
870  startvalue++;
871  break;
872  case LibUtilities::eModified_A:
873  startvalue+=2;
874  break;
875  default:
876  ASSERTL0(0,"Mapping array is not defined for this expansion");
877  break;
878  }
879 
880  for(i = 0; i < nummodes1-2; i++)
881  {
882  for(j = 0; j < nummodes0-2; j++)
883  {
884  outarray[cnt++]=startvalue+j;
885  }
886  startvalue+=nummodes0;
887  }
888  }
889 
890  int StdQuadExp::v_GetVertexMap(int localVertexId, bool useCoeffPacking)
891  {
892  int localDOF = 0;
893 
894  if(useCoeffPacking == true)
895  {
896  switch(localVertexId)
897  {
898  case 0:
899  {
900  localDOF = 0;
901  }
902  break;
903  case 1:
904  {
905  if(m_base[0]->GetBasisType()==LibUtilities::eGLL_Lagrange)
906  {
907  localDOF = m_base[0]->GetNumModes()-1;
908  }
909  else
910  {
911  localDOF = 1;
912  }
913  }
914  break;
915  case 2:
916  {
917  if(m_base[0]->GetBasisType()==LibUtilities::eGLL_Lagrange)
918  {
919  localDOF = m_base[0]->GetNumModes() * (m_base[1]->GetNumModes()-1);
920  }
921  else
922  {
923  localDOF = m_base[0]->GetNumModes();
924  }
925  }
926  break;
927  case 3:
928  {
929  if(m_base[0]->GetBasisType()==LibUtilities::eGLL_Lagrange)
930  {
931  localDOF = m_base[0]->GetNumModes()*m_base[1]->GetNumModes()-1;
932  }
933  else
934  {
935  localDOF = m_base[0]->GetNumModes()+1;
936  }
937  }
938  break;
939  default:
940  ASSERTL0(false,"eid must be between 0 and 3");
941  break;
942  }
943 
944  }
945  else
946  {
947  switch(localVertexId)
948  {
949  case 0:
950  {
951  localDOF = 0;
952  }
953  break;
954  case 1:
955  {
956  if(m_base[0]->GetBasisType()==LibUtilities::eGLL_Lagrange)
957  {
958  localDOF = m_base[0]->GetNumModes()-1;
959  }
960  else
961  {
962  localDOF = 1;
963  }
964  }
965  break;
966  case 2:
967  {
968  if(m_base[0]->GetBasisType()==LibUtilities::eGLL_Lagrange)
969  {
970  localDOF = m_base[0]->GetNumModes()*m_base[1]->GetNumModes()-1;
971  }
972  else
973  {
974  localDOF = m_base[0]->GetNumModes()+1;
975  }
976  }
977  break;
978  case 3:
979  {
980  if(m_base[0]->GetBasisType()==LibUtilities::eGLL_Lagrange)
981  {
982  localDOF = m_base[0]->GetNumModes() * (m_base[1]->GetNumModes()-1);
983  }
984  else
985  {
986  localDOF = m_base[0]->GetNumModes();
987  }
988  }
989  break;
990  default:
991  ASSERTL0(false,"eid must be between 0 and 3");
992  break;
993  }
994  }
995  return localDOF;
996  }
997 
998  void StdQuadExp::v_GetEdgeInteriorMap(const int eid,
999  const Orientation edgeOrient,
1000  Array<OneD, unsigned int> &maparray,
1001  Array<OneD, int> &signarray)
1002  {
1003  int i;
1004  const int nummodes0 = m_base[0]->GetNumModes();
1005  const int nummodes1 = m_base[1]->GetNumModes();
1006  const int nEdgeIntCoeffs = GetEdgeNcoeffs(eid)-2;
1007  const LibUtilities::BasisType bType = GetEdgeBasisType(eid);
1008 
1009  if(maparray.num_elements() != nEdgeIntCoeffs)
1010  {
1011  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1012  }
1013 
1014  if(signarray.num_elements() != nEdgeIntCoeffs)
1015  {
1016  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
1017  }
1018  else
1019  {
1020  fill( signarray.get() , signarray.get()+nEdgeIntCoeffs, 1 );
1021  }
1022 
1023  if(bType == LibUtilities::eModified_A)
1024  {
1025  switch(eid)
1026  {
1027  case 0:
1028  {
1029  for(i = 0; i < nEdgeIntCoeffs; i++)
1030  {
1031  maparray[i] = i+2;
1032  }
1033 
1034  if(edgeOrient==eBackwards)
1035  {
1036  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1037  {
1038  signarray[i] = -1;
1039  }
1040  }
1041  }
1042  break;
1043  case 1:
1044  {
1045  for(i = 0; i < nEdgeIntCoeffs; i++)
1046  {
1047  maparray[i] = (i+2)*nummodes0 + 1;
1048  }
1049 
1050  if(edgeOrient==eBackwards)
1051  {
1052  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1053  {
1054  signarray[i] = -1;
1055  }
1056  }
1057  }
1058  break;
1059  case 2:
1060  {
1061  for(i = 0; i < nEdgeIntCoeffs; i++)
1062  {
1063  maparray[i] = nummodes0+i+2;
1064  }
1065 
1066  if(edgeOrient==eForwards)
1067  {
1068  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1069  {
1070  signarray[i] = -1;
1071  }
1072  }
1073  }
1074  break;
1075  case 3:
1076  {
1077  for(i = 0; i < nEdgeIntCoeffs; i++)
1078  {
1079  maparray[i] = (i+2)*nummodes0;
1080  }
1081 
1082  if(edgeOrient==eForwards)
1083  {
1084  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1085  {
1086  signarray[i] = -1;
1087  }
1088  }
1089  }
1090  break;
1091  default:
1092  ASSERTL0(false,"eid must be between 0 and 3");
1093  break;
1094  }
1095  }
1096  else if(bType == LibUtilities::eGLL_Lagrange)
1097  {
1098  switch(eid)
1099  {
1100  case 0:
1101  {
1102  for(i = 0; i < nEdgeIntCoeffs; i++)
1103  {
1104  maparray[i] = i+1;
1105  }
1106  }
1107  break;
1108  case 1:
1109  {
1110  for(i = 0; i < nEdgeIntCoeffs; i++)
1111  {
1112  maparray[i] = (i+2)*nummodes0 - 1;
1113  }
1114  }
1115  break;
1116  case 2:
1117  {
1118  for(i = 0; i < nEdgeIntCoeffs; i++)
1119  {
1120  maparray[i] = nummodes0*nummodes1 - 2 - i;
1121  }
1122  }
1123  break;
1124  case 3:
1125  {
1126  for(i = 0; i < nEdgeIntCoeffs; i++)
1127  {
1128  maparray[i] = nummodes0*(nummodes1-2-i);
1129  }
1130  }
1131  break;
1132  default:
1133  ASSERTL0(false,"eid must be between 0 and 3");
1134  break;
1135  }
1136  if(edgeOrient == eBackwards)
1137  {
1138  reverse( maparray.get() , maparray.get()+nEdgeIntCoeffs );
1139  }
1140  }
1141  else
1142  {
1143  ASSERTL0(false,"Mapping not defined for this type of basis");
1144  }
1145 
1146  }
1147 
1148  void StdQuadExp::v_GetEdgeToElementMap(const int eid,
1149  const Orientation edgeOrient,
1150  Array<OneD, unsigned int> &maparray,
1151  Array<OneD, int> &signarray)
1152  {
1153  int i;
1154  const int nummodes0 = m_base[0]->GetNumModes();
1155  const int nummodes1 = m_base[1]->GetNumModes();
1156  const int nEdgeCoeffs = GetEdgeNcoeffs(eid);
1157  const LibUtilities::BasisType bType = GetEdgeBasisType(eid);
1158 
1159  if(maparray.num_elements() != nEdgeCoeffs)
1160  {
1161  maparray = Array<OneD, unsigned int>(nEdgeCoeffs);
1162  }
1163 
1164  if(signarray.num_elements() != nEdgeCoeffs)
1165  {
1166  signarray = Array<OneD, int>(nEdgeCoeffs,1);
1167  }
1168  else
1169  {
1170  fill( signarray.get() , signarray.get()+nEdgeCoeffs, 1 );
1171  }
1172 
1173  if(bType == LibUtilities::eModified_A)
1174  {
1175  switch(eid)
1176  {
1177  case 0:
1178  {
1179  for(i = 0; i < nEdgeCoeffs; i++)
1180  {
1181  maparray[i] = i;
1182  }
1183 
1184  if(edgeOrient==eBackwards)
1185  {
1186  swap( maparray[0] , maparray[1] );
1187 
1188  for(i = 3; i < nEdgeCoeffs; i+=2)
1189  {
1190  signarray[i] = -1;
1191  }
1192  }
1193  }
1194  break;
1195  case 1:
1196  {
1197  for(i = 0; i < nEdgeCoeffs; i++)
1198  {
1199  maparray[i] = i*nummodes0 + 1;
1200  }
1201 
1202  if(edgeOrient==eBackwards)
1203  {
1204  swap( maparray[0] , maparray[1] );
1205 
1206  for(i = 3; i < nEdgeCoeffs; i+=2)
1207  {
1208  signarray[i] = -1;
1209  }
1210  }
1211  }
1212  break;
1213  case 2:
1214  {
1215  for(i = 0; i < nEdgeCoeffs; i++)
1216  {
1217  maparray[i] = nummodes0+i;
1218  }
1219 
1220  if(edgeOrient==eForwards)
1221  {
1222  swap( maparray[0] , maparray[1] );
1223 
1224  for(i = 3; i < nEdgeCoeffs; i+=2)
1225  {
1226  signarray[i] = -1;
1227  }
1228  }
1229  }
1230  break;
1231  case 3:
1232  {
1233  for(i = 0; i < nEdgeCoeffs; i++)
1234  {
1235  maparray[i] = i*nummodes0;
1236  }
1237 
1238  if(edgeOrient==eForwards)
1239  {
1240  swap( maparray[0] , maparray[1] );
1241 
1242  for(i = 3; i < nEdgeCoeffs; i+=2)
1243  {
1244  signarray[i] = -1;
1245  }
1246  }
1247  }
1248  break;
1249  default:
1250  ASSERTL0(false,"eid must be between 0 and 3");
1251  break;
1252  }
1253  }
1254  else if(bType == LibUtilities::eGLL_Lagrange ||
1255  bType == LibUtilities::eGauss_Lagrange)
1256  {
1257  switch(eid)
1258  {
1259  case 0:
1260  {
1261  for(i = 0; i < nEdgeCoeffs; i++)
1262  {
1263  maparray[i] = i;
1264  }
1265  }
1266  break;
1267  case 1:
1268  {
1269  for(i = 0; i < nEdgeCoeffs; i++)
1270  {
1271  maparray[i] = (i+1)*nummodes0 - 1;
1272  }
1273  }
1274  break;
1275  case 2:
1276  {
1277  for(i = 0; i < nEdgeCoeffs; i++)
1278  {
1279  maparray[i] = nummodes0*nummodes1 - 1 - i;
1280  }
1281  }
1282  break;
1283  case 3:
1284  {
1285  for(i = 0; i < nEdgeCoeffs; i++)
1286  {
1287  maparray[i] = nummodes0*(nummodes1-1-i);
1288  }
1289  }
1290  break;
1291  default:
1292  ASSERTL0(false,"eid must be between 0 and 3");
1293  break;
1294  }
1295  if(edgeOrient == eBackwards)
1296  {
1297  reverse( maparray.get() , maparray.get()+nEdgeCoeffs );
1298  }
1299  }
1300  else
1301  {
1302  ASSERTL0(false,"Mapping not defined for this type of basis");
1303  }
1304  }
1305 
1306  ///////////////////////
1307  // Wrapper Functions //
1308  ///////////////////////
1309 
1310  DNekMatSharedPtr StdQuadExp::v_GenMatrix(const StdMatrixKey &mkey)
1311  {
1312  int i;
1313  int order0 = GetBasisNumModes(0);
1314  int order1 = GetBasisNumModes(1);
1315  MatrixType mtype = mkey.GetMatrixType();
1316 
1317  DNekMatSharedPtr Mat;
1318 
1319  switch(mtype)
1320  {
1321  case ePhysInterpToEquiSpaced:
1322  {
1323  int nq0 = m_base[0]->GetNumPoints();
1324  int nq1 = m_base[1]->GetNumPoints();
1325  int nq = max(nq0,nq1);
1326  int neq = LibUtilities::StdQuadData::
1327  getNumberOfCoefficients(nq, nq);
1328  Array<OneD, Array<OneD, NekDouble> > coords(neq);
1329  Array<OneD, NekDouble> coll (2);
1330  Array<OneD, DNekMatSharedPtr> I (2);
1331  Array<OneD, NekDouble> tmp (nq0);
1332 
1333  Mat = MemoryManager<DNekMat>::
1334  AllocateSharedPtr(neq, nq0 * nq1);
1335  int cnt = 0;
1336 
1337  for(int i = 0; i < nq; ++i)
1338  {
1339  for(int j = 0; j < nq; ++j,++cnt)
1340  {
1341  coords[cnt] = Array<OneD, NekDouble>(2);
1342  coords[cnt][0] = -1.0 + 2*j/(NekDouble)(nq-1);
1343  coords[cnt][1] = -1.0 + 2*i/(NekDouble)(nq-1);
1344  }
1345  }
1346 
1347  for(int i = 0; i < neq; ++i)
1348  {
1349  LocCoordToLocCollapsed(coords[i],coll);
1350 
1351  I[0] = m_base[0]->GetI(coll);
1352  I[1] = m_base[1]->GetI(coll+1);
1353 
1354  // interpolate first coordinate direction
1355  for (int j = 0; j < nq1; ++j)
1356  {
1357  NekDouble fac = (I[1]->GetPtr())[j];
1358  Vmath::Smul(nq0,fac,I[0]->GetPtr(),1,tmp,1);
1359 
1360  Vmath::Vcopy(nq0, &tmp[0], 1,
1361  Mat->GetRawPtr()+j*nq0*neq+i,neq);
1362  }
1363 
1364  }
1365  break;
1366  }
1367  case eMass:
1368  {
1369  Mat = StdExpansion::CreateGeneralMatrix(mkey);
1370  // For Fourier basis set the imaginary component of mean mode
1371  // to have a unit diagonal component in mass matrix
1372  if(m_base[0]->GetBasisType() == LibUtilities::eFourier)
1373  {
1374  for(i = 0; i < order1; ++i)
1375  {
1376  (*Mat)(order0*i+1,i*order0+1) = 1.0;
1377  }
1378  }
1379 
1380  if(m_base[1]->GetBasisType() == LibUtilities::eFourier)
1381  {
1382  for(i = 0; i < order0; ++i)
1383  {
1384  (*Mat)(order0+i ,order0+i) = 1.0;
1385  }
1386  }
1387  break;
1388  }
1389  case eFwdTrans:
1390  {
1391  Mat = MemoryManager<DNekMat>::AllocateSharedPtr(m_ncoeffs,m_ncoeffs);
1392  StdMatrixKey iprodkey(eIProductWRTBase,DetShapeType(),*this);
1393  DNekMat &Iprod = *GetStdMatrix(iprodkey);
1394  StdMatrixKey imasskey(eInvMass,DetShapeType(),*this);
1395  DNekMat &Imass = *GetStdMatrix(imasskey);
1396 
1397  (*Mat) = Imass*Iprod;
1398  break;
1399  }
1400  case eGaussDG:
1401  {
1402  ConstFactorMap factors = mkey.GetConstFactors();
1403 
1404  int edge = (int)factors[StdRegions::eFactorGaussEdge];
1405  int dir = (edge + 1) % 2;
1406  int nCoeffs = m_base[dir]->GetNumModes();
1407 
1408  const LibUtilities::PointsKey BS_p(
1409  nCoeffs, LibUtilities::eGaussGaussLegendre);
1410  const LibUtilities::BasisKey BS_k(
1411  LibUtilities::eGauss_Lagrange, nCoeffs, BS_p);
1412 
1413  Array<OneD, NekDouble> coords(1, 0.0);
1414  coords[0] = (edge == 0 || edge == 3) ? -1.0 : 1.0;
1415 
1416  LibUtilities::BasisSharedPtr basis =
1417  LibUtilities::BasisManager()[BS_k];
1418  DNekMatSharedPtr m_Ix = basis->GetI(coords);
1419 
1420  Mat = MemoryManager<DNekMat>::AllocateSharedPtr(
1421  1.0, nCoeffs);
1422  Vmath::Vcopy(nCoeffs, m_Ix->GetPtr(), 1, Mat->GetPtr(), 1);
1423  break;
1424  }
1425  default:
1426  {
1427  Mat = StdExpansion::CreateGeneralMatrix(mkey);
1428  break;
1429  }
1430  }
1431 
1432  return Mat;
1433  }
1434 
1435  DNekMatSharedPtr StdQuadExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
1436  {
1437  return GenMatrix(mkey);
1438  }
1439 
1440  ///////////////////////////////////
1441  // Operator evaluation functions //
1442  ///////////////////////////////////
1443 
1444  void StdQuadExp::v_GeneralMatrixOp_MatOp(
1445  const Array<OneD, const NekDouble> &inarray,
1446  Array<OneD,NekDouble> &outarray,
1447  const StdMatrixKey &mkey)
1448  {
1449 
1450  DNekMatSharedPtr mat = m_stdMatrixManager[mkey];
1451 
1452  if(inarray.get() == outarray.get())
1453  {
1454  Array<OneD,NekDouble> tmp(m_ncoeffs);
1455  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
1456 
1457  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, 1.0, mat->GetPtr().get(),
1458  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
1459  }
1460  else
1461  {
1462  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, 1.0, mat->GetPtr().get(),
1463  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
1464  }
1465  }
1466 
1467 
1468  void StdQuadExp::v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,
1469  const StdMatrixKey &mkey)
1470  {
1471  // Generate an orthonogal expansion
1472  int qa = m_base[0]->GetNumPoints();
1473  int qb = m_base[1]->GetNumPoints();
1474  int nmodes_a = m_base[0]->GetNumModes();
1475  int nmodes_b = m_base[1]->GetNumModes();
1476  // Declare orthogonal basis.
1477  LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
1478  LibUtilities::PointsKey pb(qb,m_base[1]->GetPointsType());
1479 
1480  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A,nmodes_a,pa);
1481  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A,nmodes_b,pb);
1482  StdQuadExp OrthoExp(Ba,Bb);
1483 
1484  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
1485 
1486  //for the "old" implementation
1487  int cutoff = (int) (mkey.GetConstFactor(eFactorSVVCutoffRatio)*min(nmodes_a,nmodes_b));
1488 
1489  //SVV parameters loaded from the .xml case file
1490  NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDiffCoeff);
1491 
1492  // project onto modal space.
1493  OrthoExp.FwdTrans(array,orthocoeffs);
1494 
1495  //counters for scanning through orthocoeffs array
1496  int j, k, cnt = 0;
1497  int nmodes = min(nmodes_a,nmodes_b);
1498 
1499  //------"New" Version August 22nd '13--------------------
1500  for(j = 0; j < nmodes_a; ++j)
1501  {
1502  for(k = 0; k < nmodes_b; ++k)
1503  {
1504  if(j + k >= cutoff) //to filter out only the "high-modes"
1505  {
1506  orthocoeffs[j*nmodes_b+k] *=
1507  (1.0+SvvDiffCoeff*exp(-(j+k-nmodes)*(j+k-nmodes)/
1508  ((NekDouble)((j+k-cutoff+1)*
1509  (j+k-cutoff+1)))));
1510  }
1511  cnt++;
1512  }
1513  }
1514 
1515  // backward transform to physical space
1516  OrthoExp.BwdTrans(orthocoeffs,array);
1517  }
1518 
1519  void StdQuadExp::v_ReduceOrderCoeffs(
1520  int numMin,
1521  const Array<OneD, const NekDouble> &inarray,
1522  Array<OneD, NekDouble> &outarray)
1523  {
1524  int n_coeffs = inarray.num_elements();
1525 
1526 
1527  Array<OneD, NekDouble> coeff(n_coeffs);
1528  Array<OneD, NekDouble> coeff_tmp(n_coeffs,0.0);
1529  Array<OneD, NekDouble> tmp;
1530  Array<OneD, NekDouble> tmp2;
1531 
1532  int nmodes0 = m_base[0]->GetNumModes();
1533  int nmodes1 = m_base[1]->GetNumModes();
1534  int numMax = nmodes0;
1535 
1536  Vmath::Vcopy(n_coeffs,inarray,1,coeff_tmp,1);
1537 
1538  const LibUtilities::PointsKey Pkey0(
1539  nmodes0, LibUtilities::eGaussLobattoLegendre);
1540  const LibUtilities::PointsKey Pkey1(
1541  nmodes1, LibUtilities::eGaussLobattoLegendre);
1542 
1543  LibUtilities::BasisKey b0(m_base[0]->GetBasisType(),nmodes0,Pkey0);
1544  LibUtilities::BasisKey b1(m_base[1]->GetBasisType(),nmodes1,Pkey1);
1545 
1546  LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A,nmodes0,Pkey0);
1547  LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_A,nmodes1,Pkey1);
1548 
1549  LibUtilities::InterpCoeff2D(
1550  b0, b1, coeff_tmp, bortho0, bortho1, coeff);
1551 
1552  Vmath::Zero(n_coeffs,coeff_tmp,1);
1553 
1554  int cnt = 0;
1555  for (int i = 0; i < numMin+1; ++i)
1556  {
1557  Vmath::Vcopy(numMin,
1558  tmp = coeff+cnt,1,
1559  tmp2 = coeff_tmp+cnt,1);
1560 
1561  cnt = i*numMax;
1562  }
1563 
1564  LibUtilities::InterpCoeff2D(
1565  bortho0, bortho1, coeff_tmp, b0, b1, outarray);
1566  }
1567 
1568 
1569  void StdQuadExp::v_MassMatrixOp(
1570  const Array<OneD, const NekDouble> &inarray,
1571  Array<OneD,NekDouble> &outarray,
1572  const StdMatrixKey &mkey)
1573  {
1574  StdExpansion::MassMatrixOp_MatFree(inarray,outarray,mkey);
1575  }
1576 
1577  void StdQuadExp::v_LaplacianMatrixOp(
1578  const Array<OneD, const NekDouble> &inarray,
1579  Array<OneD,NekDouble> &outarray,
1580  const StdMatrixKey &mkey)
1581  {
1582  StdQuadExp::v_LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
1583  }
1584 
1585  void StdQuadExp::v_LaplacianMatrixOp(const int k1, const int k2,
1586  const Array<OneD, const NekDouble> &inarray,
1587  Array<OneD,NekDouble> &outarray,
1588  const StdMatrixKey &mkey)
1589  {
1590  StdExpansion::LaplacianMatrixOp_MatFree(k1,k2,inarray,outarray,mkey);
1591  }
1592 
1593  void StdQuadExp::v_WeakDerivMatrixOp(const int i,
1594  const Array<OneD, const NekDouble> &inarray,
1595  Array<OneD,NekDouble> &outarray,
1596  const StdMatrixKey &mkey)
1597  {
1598  StdExpansion::WeakDerivMatrixOp_MatFree(i,inarray,outarray,mkey);
1599  }
1600 
1601  void StdQuadExp::v_HelmholtzMatrixOp(const Array<OneD, const NekDouble> &inarray,
1602  Array<OneD,NekDouble> &outarray,
1603  const StdMatrixKey &mkey)
1604  {
1605  StdQuadExp::v_HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
1606  }
1607 
1608  //up to here
1609  void StdQuadExp::v_MultiplyByStdQuadratureMetric(
1610  const Array<OneD, const NekDouble> &inarray,
1611  Array<OneD, NekDouble> &outarray)
1612  {
1613  int i;
1614  int nquad0 = m_base[0]->GetNumPoints();
1615  int nquad1 = m_base[1]->GetNumPoints();
1616 
1617  const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
1618  const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
1619 
1620  // multiply by integration constants
1621  for(i = 0; i < nquad1; ++i)
1622  {
1623  Vmath::Vmul(nquad0, inarray.get()+i*nquad0,1,
1624  w0.get(),1,outarray.get()+i*nquad0,1);
1625  }
1626 
1627  for(i = 0; i < nquad0; ++i)
1628  {
1629  Vmath::Vmul(nquad1,outarray.get()+i,nquad0,w1.get(),1,
1630  outarray.get()+i,nquad0);
1631  }
1632  }
1633 
1634  void StdQuadExp::v_GetSimplexEquiSpacedConnectivity(
1635  Array<OneD, int> &conn,
1636  bool standard)
1637  {
1638  int np1 = m_base[0]->GetNumPoints();
1639  int np2 = m_base[1]->GetNumPoints();
1640  int np = max(np1,np2);
1641 
1642  conn = Array<OneD, int>(6*(np-1)*(np-1));
1643 
1644  int row = 0;
1645  int rowp1 = 0;
1646  int cnt = 0;
1647  for(int i = 0; i < np-1; ++i)
1648  {
1649  rowp1 += np;
1650  for(int j = 0; j < np-1; ++j)
1651  {
1652  conn[cnt++] = row +j;
1653  conn[cnt++] = row +j+1;
1654  conn[cnt++] = rowp1 +j;
1655 
1656  conn[cnt++] = rowp1 +j+1;
1657  conn[cnt++] = rowp1 +j;
1658  conn[cnt++] = row +j+1;
1659  }
1660  row += np;
1661  }
1662  }
1663 
1664  } //end namespace
1665 }//end namespace
Generated on Wed Mar 4 2015 16:33:18 for Nektar++ by   doxygen 1.8.1.2