Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
StdSegExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File StdSegExp.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: Routines within Standard Segment Expansions
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <StdRegions/StdSegExp.h>
37 
38 namespace Nektar
39 {
40  namespace StdRegions
41  {
42 
43  /** \brief Default constructor */
44 
46  {
47  }
48 
49 
50  /** \brief Constructor using BasisKey class for quadrature points and
51  * order definition
52  *
53  * \param Ba BasisKey class definition containing order and quadrature
54  * points.
55  */
56 
58  StdExpansion(Ba.GetNumModes(), 1, Ba),
59  StdExpansion1D(Ba.GetNumModes(),Ba)
60  {
61  }
62 
63 
64  /** \brief Copy Constructor */
65 
67  StdExpansion(T),
69  {
70  }
71 
72 
74  {
75  }
76 
77  /** \brief Return Shape of region, using ShapeType enum list.
78  * i.e. Segment
79  */
81  {
83  }
84 
86  {
87 
88  bool returnval = false;
89 
91  {
92  returnval = true;
93  }
94 
96  {
97  returnval = true;
98  }
99 
100  return returnval;
101  }
102 
103 
104 
105 
106  //---------------------------------------------------------------------
107  // Integration Methods
108  //---------------------------------------------------------------------
109 
110  /** \brief Integrate the physical point list \a inarray over region
111  * and return the value
112  *
113  * \param inarray definition of function to be integrated evauluated at
114  * quadrature point of expansion.
115  * \return returns \f$\int^1_{-1} u(\xi_1)d \xi_1 \f$ where \f$inarray[i]
116  * = u(\xi_{1i}) \f$
117  */
118  NekDouble StdSegExp::v_Integral(const Array<OneD, const NekDouble>& inarray )
119  {
120  NekDouble Int = 0.0;
121  int nquad0 = m_base[0]->GetNumPoints();
122  Array<OneD, NekDouble> tmp(nquad0);
123  Array<OneD, const NekDouble> z = m_base[0]->GetZ();
124  Array<OneD, const NekDouble> w0 = m_base[0]->GetW();
125 
126  // multiply by integration constants
127  Vmath::Vmul(nquad0, inarray, 1, w0, 1, tmp, 1);
128 
129  Int = Vmath::Vsum(nquad0, tmp, 1);
130 
131  return Int;
132  }
133 
134 
135 
136 
137  //---------------------------------------------------------------------
138  // Differentiation Methods
139  //---------------------------------------------------------------------
140 
141 
142  /** \brief Evaluate the derivative \f$ d/d{\xi_1} \f$ at the physical
143  * quadrature points given by \a inarray and return in \a outarray.
144  *
145  * This is a wrapper around StdExpansion1D::Tensor_Deriv
146  * \param inarray array of a function evaluated at the quadrature points
147  * \param outarray the resulting array of the derivative \f$
148  * du/d_{\xi_1}|_{\xi_{1i}} \f$ will be stored in the array \a outarra
149  */
150 
151  void StdSegExp::v_PhysDeriv(const Array<OneD, const NekDouble>& inarray,
152  Array<OneD, NekDouble> &out_d0,
153  Array<OneD, NekDouble> &out_d1,
154  Array<OneD, NekDouble> &out_d2)
155  {
156  PhysTensorDeriv(inarray,out_d0);
157  }
158 
159 
160  void StdSegExp::v_PhysDeriv(const int dir,
161  const Array<OneD, const NekDouble>& inarray,
162  Array<OneD, NekDouble> &outarray)
163  {
164  ASSERTL1(dir==0,"input dir is out of range");
165  PhysTensorDeriv(inarray,outarray);
166  // PhysDeriv(inarray, outarray);
167  }
168 
170  const Array<OneD, const NekDouble>& inarray,
171  Array<OneD, NekDouble> &out_d0,
172  Array<OneD, NekDouble> &out_d1,
173  Array<OneD, NekDouble> &out_d2)
174  {
175  PhysTensorDeriv(inarray,out_d0);
176  // PhysDeriv(inarray, out_d0);
177  }
178 
180  const int dir,
181  const Array<OneD, const NekDouble>& inarray,
182  Array<OneD, NekDouble> &outarray)
183  {
184  ASSERTL1(dir==0,"input dir is out of range");
185  PhysTensorDeriv(inarray,outarray);
186  // PhysDeriv(inarray, outarray);
187  }
188 
189 
190 
191  //---------------------------------------------------------------------
192  // Transforms
193  //---------------------------------------------------------------------
194 
195  /** \brief Backward transform from coefficient space given
196  * in \a inarray and evaluate at the physical quadrature
197  * points \a outarray
198  *
199  * Operation can be evaluated as \f$ u(\xi_{1i}) =
200  * \sum_{p=0}^{order-1} \hat{u}_p \phi_p(\xi_{1i}) \f$ or equivalently
201  * \f$ {\bf u} = {\bf B}^T {\bf \hat{u}} \f$ where
202  * \f${\bf B}[i][j] = \phi_i(\xi_{1j}), \mbox{\_coeffs}[p] = {\bf
203  * \hat{u}}[p] \f$
204  *
205  * The function takes the coefficient array \a inarray as
206  * input for the transformation
207  *
208  * \param inarray: the coeffficients of the expansion
209  *
210  * \param outarray: the resulting array of the values of the function at
211  * the physical quadrature points will be stored in the array \a outarray
212  */
213 
215  const Array<OneD, const NekDouble>& inarray,
216  Array<OneD, NekDouble> &outarray)
217  {
218  int nquad = m_base[0]->GetNumPoints();
219 
220  if(m_base[0]->Collocation())
221  {
222  Vmath::Vcopy(nquad, inarray, 1, outarray, 1);
223  }
224  else
225  {
226 
227 #ifdef NEKTAR_USING_DIRECT_BLAS_CALLS
228 
229  Blas::Dgemv('N',nquad,m_base[0]->GetNumModes(),1.0, (m_base[0]->GetBdata()).get(),
230  nquad,&inarray[0],1,0.0,&outarray[0],1);
231 
232 #else //NEKTAR_USING_DIRECT_BLAS_CALLS
233 
235  NekVector<NekDouble> out(nquad,outarray,eWrapper);
236  NekMatrix<NekDouble> B(nquad,m_ncoeffs,m_base[0]->GetBdata(),eWrapper);
237  out = B * in;
238 
239 #endif //NEKTAR_USING_DIRECT_BLAS_CALLS
240 
241  }
242  }
243 
244  /** \brief Forward transform from physical quadrature space stored in
245  * \a inarray and evaluate the expansion coefficients and store in
246  * \a outarray
247  *
248  * Perform a forward transform using a Galerkin projection by taking the
249  * inner product of the physical points and multiplying by the inverse of
250  * the mass matrix using the Solve method of the standard matrix
251  * container holding the local mass matrix, i.e. \f$ {\bf \hat{u}} =
252  * {\bf M}^{-1} {\bf I} \f$ where \f$ {\bf I}[p] = \int^1_{-1}
253  * \phi_p(\xi_1) u(\xi_1) d\xi_1 \f$
254  *
255  * This function stores the expansion coefficients calculated by the
256  * transformation in the coefficient space array \a outarray
257  *
258  * \param inarray: array of physical quadrature points to be transformed
259  *
260  * \param outarray: the coeffficients of the expansion
261  */
262 
263  void StdSegExp::v_FwdTrans(const Array<OneD, const NekDouble>& inarray,
264  Array<OneD, NekDouble> &outarray)
265  {
266  if(m_base[0]->Collocation())
267  {
268  Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
269  }
270  else
271  {
272  v_IProductWRTBase(inarray,outarray);
273 
274  // get Mass matrix inverse
275  StdMatrixKey masskey(eInvMass,v_DetShapeType(),*this);
276  DNekMatSharedPtr matsys = GetStdMatrix(masskey);
277 
278  NekVector<NekDouble> in(m_ncoeffs,outarray,eCopy);
280 
281  out = (*matsys)*in;
282  }
283  }
284 
286  const Array<OneD, const NekDouble>& inarray,
287  Array<OneD, NekDouble> &outarray)
288  {
289  if(m_base[0]->Collocation())
290  {
291  Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
292  }
293  else
294  {
295  int nInteriorDofs = m_ncoeffs-2;
296  int offset;
297 
298  switch(m_base[0]->GetBasisType())
299  {
301  {
302  offset = 1;
303  }
304  break;
306  {
307  nInteriorDofs = m_ncoeffs;
308  offset = 0;
309  }
310  break;
313  {
314  offset = 2;
315  }
316  break;
317  default:
318  ASSERTL0(false,"This type of FwdTrans is not defined for this expansion type");
319  }
320 
321  fill(outarray.get(), outarray.get()+m_ncoeffs, 0.0 );
322 
324  {
325  outarray[GetVertexMap(0)] = inarray[0];
326  outarray[GetVertexMap(1)] = inarray[m_base[0]->GetNumPoints()-1];
327 
328  if(m_ncoeffs>2)
329  {
330  // ideally, we would like to have tmp0 to be replaced by
331  // outarray (currently MassMatrixOp does not allow aliasing)
332  Array<OneD, NekDouble> tmp0(m_ncoeffs);
333  Array<OneD, NekDouble> tmp1(m_ncoeffs);
334 
335  StdMatrixKey masskey(eMass,v_DetShapeType(),*this);
336  MassMatrixOp(outarray,tmp0,masskey);
337  v_IProductWRTBase(inarray,tmp1);
338 
339  Vmath::Vsub(m_ncoeffs, tmp1, 1, tmp0, 1, tmp1, 1);
340 
341  // get Mass matrix inverse (only of interior DOF)
342  DNekMatSharedPtr matsys =
343  (m_stdStaticCondMatrixManager[masskey])-> GetBlock(1,1);
344 
345  Blas::Dgemv('N',nInteriorDofs,nInteriorDofs,1.0,
346  &(matsys->GetPtr())[0],nInteriorDofs,tmp1.get()
347  +offset,1,0.0,outarray.get()+offset,1);
348  }
349  }
350  else
351  {
352  StdSegExp::v_FwdTrans(inarray, outarray);
353  }
354  }
355 
356  }
357 
358 
359  void StdSegExp::v_BwdTrans_SumFac(const Array<OneD, const NekDouble>& inarray,
360  Array<OneD, NekDouble> &outarray)
361  {
362  v_BwdTrans(inarray, outarray);
363  }
364 
365 
366 
367  //---------------------------------------------------------------------
368  // Inner product functions
369  //---------------------------------------------------------------------
370 
371 
372 
373  /** \brief Inner product of \a inarray over region with respect to
374  * expansion basis \a base and return in \a outarray
375  *
376  * Calculate \f$ I[p] = \int^{1}_{-1} \phi_p(\xi_1) u(\xi_1) d\xi_1
377  * = \sum_{i=0}^{nq-1} \phi_p(\xi_{1i}) u(\xi_{1i}) w_i \f$ where
378  * \f$ outarray[p] = I[p], inarray[i] = u(\xi_{1i}), base[p*nq+i] =
379  * \phi_p(\xi_{1i}) \f$.
380  *
381  * \param base an array defining the local basis for the inner product
382  * usually passed from Basis->GetBdata() or Basis->GetDbdata()
383  * \param inarray: physical point array of function to be integrated
384  * \f$ u(\xi_1) \f$
385  * \param coll_check flag to identify when a Basis->Collocation() call
386  * should be performed to see if this is a GLL_Lagrange basis with a
387  * collocation property. (should be set to 0 if taking the inner
388  * product with respect to the derivative of basis)
389  * \param outarray the values of the inner product with respect to
390  * each basis over region will be stored in the array \a outarray as
391  * output of the function
392  */
393 
395  const Array<OneD, const NekDouble> &base,
396  const Array<OneD, const NekDouble> &inarray,
397  Array<OneD, NekDouble> &outarray,
398  int coll_check)
399  {
400  int nquad = m_base[0]->GetNumPoints();
401  Array<OneD, NekDouble> tmp(nquad);
402  Array<OneD, const NekDouble> z = m_base[0]->GetZ();
403  Array<OneD, const NekDouble> w = m_base[0]->GetW();
404 
405  Vmath::Vmul(nquad, inarray, 1, w, 1, tmp, 1);
406 
407  /* Comment below was a bug for collocated basis
408  if(coll_check&&m_base[0]->Collocation())
409  {
410  Vmath::Vcopy(nquad, tmp, 1, outarray, 1);
411  }
412  else
413  {
414  Blas::Dgemv('T',nquad,m_ncoeffs,1.0,base.get(),nquad,
415  &tmp[0],1,0.0,outarray.get(),1);
416  }*/
417 
418  // Correct implementation
419  Blas::Dgemv('T',nquad,m_ncoeffs,1.0,base.get(),nquad,
420  &tmp[0],1,0.0,outarray.get(),1);
421  }
422 
423  /** \brief Inner product of \a inarray over region with respect to the
424  * expansion basis (this)->m_base[0] and return in \a outarray
425  *
426  * Wrapper call to StdSegExp::IProductWRTBase
427  * \param inarray array of function values evaluated at the physical
428  * collocation points
429  * \param outarray the values of the inner product with respect to
430  * each basis over region will be stored in the array \a outarray as
431  * output of the function
432  */
433  void StdSegExp::v_IProductWRTBase(const Array<OneD, const NekDouble>& inarray,
434  Array<OneD, NekDouble> &outarray)
435  {
436  v_IProductWRTBase(m_base[0]->GetBdata(),inarray,outarray,1);
437  }
438 
440  const int dir,
441  const Array<OneD, const NekDouble>& inarray,
442  Array<OneD, NekDouble> & outarray)
443  {
444  ASSERTL1(dir >= 0 && dir < 1,"input dir is out of range");
445  v_IProductWRTBase(m_base[0]->GetDbdata(),inarray,outarray,1);
446  }
447 
449  const Array<OneD, const NekDouble>& inarray,
450  Array<OneD, NekDouble> &outarray)
451  {
452  v_IProductWRTBase(m_base[0]->GetBdata(),inarray,outarray,1);
453  }
454 
455 
456 
457 
458 
459  //----------------------------
460  // Evaluation
461  //----------------------------
462 
463  void StdSegExp::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
464  {
465  int nquad = m_base[0]->GetNumPoints();
466  const NekDouble * base = m_base[0]->GetBdata().get();
467 
468  ASSERTL2(mode <= m_ncoeffs,
469  "calling argument mode is larger than total expansion order");
470 
471  Vmath::Vcopy(nquad,(NekDouble *)base+mode*nquad,1, &outarray[0],1);
472  }
473 
474 
476  const Array<OneD, const NekDouble>& coords,
477  const Array<OneD, const NekDouble>& physvals)
478  {
479  return StdExpansion1D::v_PhysEvaluate(coords, physvals);
480  }
481 
483  const Array<OneD, const NekDouble> &inarray,
484  Array<OneD, NekDouble> &outarray,
485  const StdMatrixKey &mkey)
486  {
487  int nquad = m_base[0]->GetNumPoints();
488 
489  Array<OneD,NekDouble> physValues(nquad);
490  Array<OneD,NekDouble> dPhysValuesdx(nquad);
491 
492  v_BwdTrans(inarray,physValues);
493 
494  // Laplacian matrix operation
495  v_PhysDeriv(physValues,dPhysValuesdx);
496  v_IProductWRTBase(m_base[0]->GetDbdata(),dPhysValuesdx,outarray,1);
497  }
498 
499 
501  const Array<OneD, const NekDouble> &inarray,
502  Array<OneD, NekDouble> &outarray,
503  const StdMatrixKey &mkey)
504  {
505  int nquad = m_base[0]->GetNumPoints();
506 
507  Array<OneD,NekDouble> physValues(nquad);
508  Array<OneD,NekDouble> dPhysValuesdx(nquad);
509  Array<OneD,NekDouble> wsp(m_ncoeffs);
510 
511  v_BwdTrans(inarray,physValues);
512 
513  // mass matrix operation
514  v_IProductWRTBase((m_base[0]->GetBdata()),physValues,wsp,1);
515 
516  // Laplacian matrix operation
517  v_PhysDeriv(physValues,dPhysValuesdx);
518  v_IProductWRTBase(m_base[0]->GetDbdata(),dPhysValuesdx,outarray,1);
519  Blas::Daxpy(m_ncoeffs, mkey.GetConstFactor(eFactorLambda), wsp.get(), 1, outarray.get(), 1);
520  }
521 
522 
524  Array<OneD, NekDouble> &coords_0,
525  Array<OneD, NekDouble> &coords_1,
526  Array<OneD, NekDouble> &coords_2)
527  {
528  Blas::Dcopy(GetNumPoints(0),(m_base[0]->GetZ()).get(),
529  1,&coords_0[0],1);
530  }
531 
532 
533 
534 
535 
536  //---------------------------------------------------------------------
537  // Helper functions
538  //---------------------------------------------------------------------
539 
541  {
542  return 2;
543  }
544 
546  {
547  return 2;
548  }
549 
551  {
552  return 2;
553  }
554 
556  const std::vector<unsigned int> &nummodes,
557  int &modes_offset)
558  {
559  int nmodes = nummodes[modes_offset];
560  modes_offset += 1;
561 
562  return nmodes;
563  }
564 
565  //---------------------------------------------------------------------
566  // Wrapper functions
567  //---------------------------------------------------------------------
568 
570  {
571  DNekMatSharedPtr Mat;
572  MatrixType mattype;
573 
574  switch(mattype = mkey.GetMatrixType())
575  {
576  case eFwdTrans:
577  {
579  StdMatrixKey iprodkey(eIProductWRTBase,v_DetShapeType(),*this);
580  DNekMat &Iprod = *GetStdMatrix(iprodkey);
581  StdMatrixKey imasskey(eInvMass,v_DetShapeType(),*this);
582  DNekMat &Imass = *GetStdMatrix(imasskey);
583 
584  (*Mat) = Imass*Iprod;
585  }
586  break;
587  default:
588  {
590 
591  if(mattype == eMass)
592  {
593  // For Fourier basis set the imaginary component
594  // of mean mode to have a unit diagonal component
595  // in mass matrix
597  {
598  (*Mat)(1,1) = 1.0;
599  }
600  }
601  }
602  break;
603  }
604 
605  return Mat;
606  }
607 
608 
610  {
611  return v_GenMatrix(mkey);
612  }
613 
614  //---------------------------------------------------------------------
615  // Mappings
616  //---------------------------------------------------------------------
617 
618 
619  void StdSegExp::v_GetBoundaryMap(Array<OneD, unsigned int>& outarray)
620  {
621  if(outarray.num_elements() != NumBndryCoeffs())
622  {
623  outarray = Array<OneD, unsigned int>(NumBndryCoeffs());
624  }
625  const LibUtilities::BasisType Btype = GetBasisType(0);
626  int nummodes = m_base[0]->GetNumModes();
627 
628  outarray[0] = 0;
629 
630  switch(Btype)
631  {
636  outarray[1]= nummodes-1;
637  break;
640  outarray[1] = 1;
641  break;
642  default:
643  ASSERTL0(0,"Mapping array is not defined for this expansion");
644  break;
645  }
646  }
647 
648  void StdSegExp::v_GetInteriorMap(Array<OneD, unsigned int>& outarray)
649  {
650  int i;
651  if(outarray.num_elements()!=GetNcoeffs()-NumBndryCoeffs())
652  {
653  outarray = Array<OneD, unsigned int>(GetNcoeffs()-NumBndryCoeffs());
654  }
655  const LibUtilities::BasisType Btype = GetBasisType(0);
656 
657  switch(Btype)
658  {
663  for(i = 0 ; i < GetNcoeffs()-2;i++)
664  {
665  outarray[i] = i+1;
666  }
667  break;
670  for(i = 0 ; i < GetNcoeffs()-2;i++)
671  {
672  outarray[i] = i+2;
673  }
674  break;
675  default:
676  ASSERTL0(0,"Mapping array is not defined for this expansion");
677  break;
678  }
679  }
680 
681  int StdSegExp::v_GetVertexMap(int localVertexId,bool useCoeffPacking)
682  {
683  ASSERTL0((localVertexId==0)||(localVertexId==1),"local vertex id"
684  "must be between 0 or 1");
685 
686  int localDOF = localVertexId;
687 
689  (localVertexId==1) )
690  {
691  localDOF = m_base[0]->GetNumModes()-1;
692  }
693  return localDOF;
694  }
695 
696 
697  }//end namespace
698 }//end namespace
699