Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Basis.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File StdBasis.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: Basis definition
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <LibUtilities/Foundations/Basis.h>
37 #include <LibUtilities/Foundations/ManagerAccess.h>
38 #include <LibUtilities/Foundations/Points.h>
39 #include <LibUtilities/Polylib/Polylib.h>
40 #include <LibUtilities/BasicUtils/ErrorUtil.hpp>
41 #include <LibUtilities/LinearAlgebra/Blas.hpp>
42 #include <boost/math/special_functions/gamma.hpp>
43 
44 namespace Nektar
45 {
46  namespace LibUtilities
47  {
48 
49  bool operator<(const BasisKey &lhs, const BasisKey &rhs)
50  {
51  PointsKey lhsPointsKey = lhs.GetPointsKey();
52  PointsKey rhsPointsKey = rhs.GetPointsKey();
53 
54  if (lhsPointsKey < rhsPointsKey)
55  {
56  return true;
57  }
58  if (lhsPointsKey != rhsPointsKey)
59  {
60  return false;
61  }
62 
63  if (lhs.m_nummodes < rhs.m_nummodes)
64  {
65  return true;
66  }
67  if (lhs.m_nummodes > rhs.m_nummodes)
68  {
69  return false;
70  }
71 
72  return (lhs.m_basistype < rhs.m_basistype);
73  }
74 
75  bool operator>(const BasisKey &lhs, const BasisKey &rhs)
76  {
77  return (rhs < lhs);
78  }
79 
80  bool BasisKey::opLess::operator()(const BasisKey &lhs, const BasisKey &rhs) const
81  {
82  return (lhs.m_basistype < rhs.m_basistype);
83  }
84 
85  std::ostream& operator<<(std::ostream& os, const BasisKey& rhs)
86  {
87  os << "NumModes: " << rhs.GetNumModes() << " BasisType: " << BasisTypeMap[rhs.GetBasisType()];
88  os << " " << rhs.GetPointsKey() << std::endl;
89 
90  return os;
91  }
92 
93  Basis::Basis(const BasisKey &bkey):
94  m_basisKey(bkey),
95  m_points(PointsManager()[bkey.GetPointsKey()]),
96  m_bdata(bkey.GetTotNumModes()*bkey.GetTotNumPoints()),
97  m_dbdata(bkey.GetTotNumModes()*bkey.GetTotNumPoints())
98  {
99  m_InterpManager.RegisterGlobalCreator(boost::bind(&Basis::CalculateInterpMatrix,this,_1));
100  }
101 
102  boost::shared_ptr<Basis> Basis::Create(const BasisKey &bkey)
103  {
104  boost::shared_ptr<Basis> returnval(new Basis(bkey));
105  returnval->Initialize();
106 
107  return returnval;
108  }
109 
110  void Basis::Initialize()
111 
112  {
113  ASSERTL0(GetNumModes()>0, "Cannot call Basis initialisation with zero or negative order");
114  ASSERTL0(GetTotNumPoints()>0, "Cannot call Basis initialisation with zero or negative numbers of points");
115 
116  GenBasis();
117  };
118 
119  /** \brief Calculate the interpolation Matrix for coefficient from
120  * one base (m_basisKey) to another (tbasis0)
121  */
122  boost::shared_ptr< NekMatrix<NekDouble> > Basis::CalculateInterpMatrix(const BasisKey &tbasis0)
123  {
124  int dim = m_basisKey.GetNumModes();
125  const PointsKey pkey(dim,LibUtilities::eGaussLobattoLegendre);
126  BasisKey fbkey(m_basisKey.GetBasisType(),dim,pkey);
127  BasisKey tbkey(tbasis0.GetBasisType(),dim,pkey);
128 
129  // "Constructur" of the basis
130  BasisSharedPtr fbasis = BasisManager()[fbkey];
131  BasisSharedPtr tbasis = BasisManager()[tbkey];
132 
133  // Get B Matrices
134  Array<OneD, NekDouble> fB_data = fbasis->GetBdata();
135  Array<OneD, NekDouble> tB_data = tbasis->GetBdata();
136 
137  // Convert to a NekMatrix
138  NekMatrix<NekDouble> fB(dim,dim,fB_data);
139  NekMatrix<NekDouble> tB(dim,dim,tB_data);
140 
141  // Invert the "to" matrix: tu = tB^(-1)*fB fu = ftB fu
142  tB.Invert();
143 
144  // Compute transformation matrix
145  Array<OneD, NekDouble> zero1D(dim*dim,0.0);
146  boost::shared_ptr< NekMatrix<NekDouble> > ftB(MemoryManager<NekMatrix<NekDouble> >::AllocateSharedPtr(dim,dim,zero1D));
147  (*ftB) = tB*fB;
148 
149  return ftB;
150  }
151 
152  // Method used to generate appropriate basis
153  /** The following expansions are generated depending on the
154  * enum type defined in \a m_basisKey.m_basistype:
155  *
156  * NOTE: This definition does not follow the order in the
157  * Karniadakis \& Sherwin book since this leads to a more
158  * compact hierarchical pattern for implementation
159  * purposes. The order of these modes dictates the
160  * ordering of the expansion coefficients.
161  *
162  * In the following m_numModes = P
163  *
164  * \a eModified_A:
165  *
166  * m_bdata[i + j*m_numpoints] =
167  * \f$ \phi^a_i(z_j) = \left \{
168  * \begin{array}{ll} \left ( \frac{1-z_j}{2}\right ) & i = 0 \\
169  * \\
170  * \left ( \frac{1+z_j}{2}\right ) & i = 1 \\
171  * \\
172  * \left ( \frac{1-z_j}{2}\right )\left ( \frac{1+z_j}{2}\right )
173  * P^{1,1}_{i-2}(z_j) & 2\leq i < P\\
174  * \end{array} \right . \f$
175  *
176  * \a eModified_B:
177  *
178  * m_bdata[n(i,j) + k*m_numpoints] =
179  * \f$ \phi^b_{ij}(z_k) = \left \{ \begin{array}{lll}
180  * \phi^a_j(z_k) & i = 0, & 0\leq j < P \\
181  * \\
182  * \left ( \frac{1-z_k}{2}\right )^{i} & 1 \leq i < P,& j = 0 \\
183  * \\
184  * \left ( \frac{1-z_k}{2}\right )^{i} \left ( \frac{1+z_k}{2}\right )
185  * P^{2i-1,1}_{j-1}(z_k) & 1 \leq i < P,\ & 1\leq j < P-i\ \\
186  * \end{array} \right . , \f$
187  *
188  * where \f$ n(i,j) \f$ is a consecutive ordering of the
189  * triangular indices \f$ 0 \leq i, i+j < P \f$ where \a j
190  * runs fastest.
191  *
192  *
193  * \a eModified_C:
194  *
195  * m_bdata[n(i,j,k) + l*m_numpoints] =
196  * \f$ \phi^c_{ij,k}(z_l) = \phi^b_{i+j,k}(z_l) =
197  * \left \{ \begin{array}{llll}
198  * \phi^b_{j,k}(z_l) & i = 0, & 0\leq j < P & 0\leq k < P-j\\
199  * \\
200  * \left ( \frac{1-z_l}{2}\right )^{i+j} & 1\leq i < P,\
201  * & 0\leq j < P-i,\ & k = 0 \\
202  * \\
203  * \left ( \frac{1-z_l}{2}\right )^{i+j}
204  * \left ( \frac{1+z_l}{2}\right )
205  * P^{2i+2j-1,1}_{k-1}(z_k) & 1\leq i < P,& 0\leq j < P-i&
206  * 1\leq k < P-i-j \\
207  * \\
208  * \end{array} \right . , \f$
209  *
210  * where \f$ n(i,j,k) \f$ is a consecutive ordering of the
211  * triangular indices \f$ 0 \leq i, i+j, i+j+k < P \f$ where \a k
212  * runs fastest, then \a j and finally \a i.
213  *
214  */
215  void Basis::GenBasis()
216  {
217  int i,p,q;
218  NekDouble scal;
219  Array<OneD, NekDouble> modeSharedArray;
220  NekDouble *mode;
221  Array<OneD, const NekDouble> z;
222  Array<OneD, const NekDouble> w;
223  const NekDouble *D;
224 
225  m_points->GetZW(z,w);
226 
227  D = &(m_points->GetD()->GetPtr())[0];
228  int numModes = GetNumModes();
229  int numPoints = GetNumPoints();
230 
231  switch(GetBasisType())
232  {
233 
234  /** \brief Orthogonal basis A
235 
236  \f$\tilde \psi_p^a (\eta_1) = L_p(\eta_1) = P_p^{0,0}(\eta_1)\f$
237 
238  */
239  case eOrtho_A:
240  case eLegendre:
241  mode = m_bdata.data();
242 
243  for (p=0; p<numModes; ++p, mode += numPoints)
244  {
245  Polylib::jacobfd(numPoints, z.data(), mode, NULL, p, 0.0, 0.0);
246  // normalise
247  scal = sqrt(0.5*(2.0*p+1.0));
248  for(i = 0; i < numPoints; ++i)
249  {
250  mode[i] *= scal;
251  }
252  }
253  // define derivative basis
254  Blas::Dgemm('n','n',numPoints,numModes,numPoints,1.0,D,numPoints,
255  m_bdata.data(),numPoints,0.0,m_dbdata.data(),numPoints);
256  break;
257 
258  /** \brief Orthogonal basis B
259 
260  \f$\tilde \psi_{pq}^b(\eta_2) = \left ( {1 - \eta_2} \over 2 \right)^p P_q^{2p+1,0}(\eta_2)\f$ \\
261 
262  */
263 
264  // This is tilde psi_pq in Spen's book, page 105
265  // The 3-dimensional array is laid out in memory such that
266  // 1) Eta_y values are the changing the fastest, then q and p.
267  // 2) q index increases by the stride of numPoints.
268  case eOrtho_B:
269  {
270  NekDouble *mode = m_bdata.data();
271 
272  for( int p = 0; p < numModes; ++p )
273  {
274  for( int q = 0; q < numModes - p; ++q, mode += numPoints )
275  {
276  Polylib::jacobfd(numPoints, z.data(), mode, NULL, q, 2*p + 1.0, 0.0);
277  for( int j = 0; j < numPoints; ++j )
278  {
279  mode[j] *= sqrt(p+q+1.0)*pow(0.5*(1.0 - z[j]), p);
280  }
281  }
282  }
283 
284  // define derivative basis
285  Blas::Dgemm('n','n',numPoints,numModes*(numModes+1)/2,numPoints,1.0,D,numPoints,
286  m_bdata.data(),numPoints,0.0,m_dbdata.data(),numPoints);
287  }
288  break;
289 
290  /** \brief Orthogonal basis C
291 
292  \f$\tilde \psi_{pqr}^c = \left ( {1 - \eta_3} \over 2 \right)^{p+q} P_r^{2p+2q+2, 0}(\eta_3)\f$ \\
293 
294  */
295 
296  // This is tilde psi_pqr in Spen's book, page 105
297  // The 4-dimensional array is laid out in memory such that
298  // 1) Eta_z values are the changing the fastest, then r, q, and finally p.
299  // 2) r index increases by the stride of numPoints.
300  case eOrtho_C:
301  {
302  int P = numModes - 1, Q = numModes - 1, R = numModes - 1;
303  NekDouble *mode = m_bdata.data();
304 
305  for( int p = 0; p <= P; ++p )
306  {
307  for( int q = 0; q <= Q - p; ++q )
308  {
309  for( int r = 0; r <= R - p - q; ++r, mode += numPoints )
310  {
311  Polylib::jacobfd(numPoints, z.data(), mode, NULL, r, 2*p + 2*q + 2.0, 0.0);
312  for( int k = 0; k < numPoints; ++k )
313  {
314  // Note factor of 0.5 is part of normalisation
315  mode[k] *= pow(0.5*(1.0 - z[k]), p+q);
316 
317  // finish normalisation
318  mode[k] *= sqrt(r+p+q+1.5);
319  }
320  }
321  }
322  }
323 
324  // Define derivative basis
325  Blas::Dgemm('n','n',numPoints,numModes*(numModes+1)*
326  (numModes+2)/6,numPoints,1.0, D, numPoints,
327  m_bdata.data(),numPoints,0.0,m_dbdata.data(),numPoints);
328  }
329  break;
330 
331  case eModified_A:
332 
333  // Note the following packing deviates from the
334  // definition in the Book by Karniadakis in that we
335  // put the vertex degrees of freedom at the lower
336  // index range to follow a more hierarchic structure.
337 
338  for(i = 0; i < numPoints; ++i)
339  {
340  m_bdata[i] = 0.5*(1-z[i]);
341  m_bdata[numPoints + i] = 0.5*(1+z[i]);
342  }
343 
344  mode = m_bdata.data() + 2*numPoints;
345 
346  for(p = 2; p < numModes; ++p, mode += numPoints)
347  {
348  Polylib::jacobfd(numPoints, z.data(), mode, NULL, p-2,1.0,1.0);
349 
350  for(i = 0; i < numPoints; ++i)
351  {
352  mode[i] *= m_bdata[i]*m_bdata[numPoints+i];
353  }
354  }
355 
356  // define derivative basis
357  Blas::Dgemm('n','n',numPoints,numModes,numPoints,1.0,D,
358  numPoints,m_bdata.data(),numPoints,0.0,m_dbdata.data(),
359  numPoints);
360  break;
361 
362  case eModified_B:
363  {
364 
365  // Note the following packing deviates from the
366  // definition in the Book by Karniadakis in two
367  // ways. 1) We put the vertex degrees of freedom
368  // at the lower index range to follow a more
369  // hierarchic structure. 2) We do not duplicate
370  // the singular vertex definition so that only a
371  // triangular number (i.e. (modes)*(modes+1)/2) of
372  // modes are required consistent with the
373  // orthogonal basis.
374 
375  // In the current structure the q index runs
376  // faster than the p index so that the matrix has
377  // a more compact structure
378 
379  const NekDouble *one_m_z_pow, *one_p_z;
380 
381  // bdata should be of size order*(order+1)/2*zorder
382 
383  // first fow
384  for(i = 0; i < numPoints; ++i)
385  {
386  m_bdata[0*numPoints + i] = 0.5*(1-z[i]);
387  m_bdata[1*numPoints + i] = 0.5*(1+z[i]);
388  }
389 
390  mode = m_bdata.data() + 2*numPoints;
391 
392  for(q = 2; q < numModes; ++q, mode+=numPoints)
393  {
394  Polylib::jacobfd(numPoints, z.data(), mode, NULL, q-2,1.0,1.0);
395 
396  for(i = 0; i < numPoints; ++i)
397  {
398  mode[i] *= m_bdata[i]*m_bdata[numPoints+i];
399  }
400  }
401 
402  // second row
403  for(i = 0; i < numPoints; ++i)
404  {
405  mode[i] = 0.5*(1-z[i]);
406  }
407 
408  mode += numPoints;
409 
410  for(q = 2; q < numModes; ++q, mode+=numPoints)
411  {
412  Polylib::jacobfd(numPoints, z.data(), mode, NULL, q-2,1.0,1.0);
413 
414  for(i = 0; i < numPoints; ++i)
415  {
416  mode[i] *= m_bdata[i]*m_bdata[numPoints+i];
417  }
418  }
419 
420  // third and higher rows
421  one_m_z_pow = m_bdata.data();
422  one_p_z = m_bdata.data()+numPoints;
423 
424  for(p = 2; p < numModes; ++p)
425  {
426  for(i = 0; i < numPoints; ++i)
427  {
428  mode[i] = m_bdata[i]*one_m_z_pow[i];
429  }
430 
431  one_m_z_pow = mode;
432  mode += numPoints;
433 
434  for(q = 1; q < numModes-p; ++q, mode+=numPoints)
435  {
436  Polylib::jacobfd(numPoints,z.data(),mode,NULL,q-1,2*p-1,1.0);
437 
438  for(i = 0; i < numPoints; ++i)
439  {
440  mode[i] *= one_m_z_pow[i]*one_p_z[i];
441  }
442  }
443  }
444 
445  Blas::Dgemm('n','n',numPoints,numModes*(numModes+1)/2,
446  numPoints,1.0,D,numPoints,
447  m_bdata.data(),numPoints,0.0,m_dbdata.data(),numPoints);
448  }
449  break;
450 
451 
452  case eModified_C:
453  {
454  // Note the following packing deviates from the
455  // definition in the Book by Karniadakis in two
456  // ways. 1) We put the vertex degrees of freedom
457  // at the lower index range to follow a more
458  // hierarchic structure. 2) We do not duplicate
459  // the singular vertex definition (or the
460  // duplicated face information in the book ) so
461  // that only a tetrahedral number
462  // (i.e. (modes)*(modes+1)*(modes+2)/6) of modes
463  // are required consistent with the orthogonal
464  // basis.
465 
466  // In the current structure the r index runs
467  // fastest rollowed by q and than the p index so
468  // that the matrix has a more compact structure
469 
470  // Note that eModified_C is a re-organisation/
471  // duplication of eModified_B so will get a
472  // temporary Modified_B expansion and copy the
473  // correct components.
474 
475  // Generate Modified_B basis;
476  BasisKey ModBKey(eModified_B,m_basisKey.GetNumModes(),
477  m_basisKey.GetPointsKey());
478  BasisSharedPtr ModB = BasisManager()[ModBKey];
479 
480  Array<OneD, const NekDouble> ModB_data = ModB->GetBdata();
481 
482  // Copy Modified_B basis into first
483  // (numModes*(numModes+1)/2)*numPoints entires of
484  // bdata. This fills in the complete (r,p) face.
485 
486  // Set up \phi^c_{p,q,r} = \phi^b_{p+q,r}
487 
488  int N;
489  int B_offset = 0;
490  int offset = 0;
491  for(p = 0; p < numModes; ++p)
492  {
493  N = numPoints*(numModes-p)*(numModes-p+1)/2;
494  Vmath::Vcopy(N, &ModB_data[0]+B_offset,1,&m_bdata[0] + offset,1);
495  B_offset += numPoints*(numModes-p);
496  offset += N;
497  }
498 
499  // set up derivative of basis.
500  Blas::Dgemm('n','n',numPoints,
501  numModes*(numModes+1)*(numModes+2)/6,
502  numPoints,1.0,D,numPoints,
503  m_bdata.data(),numPoints,0.0,
504  m_dbdata.data(),numPoints);
505  }
506  break;
507 
508  case eGLL_Lagrange:
509  {
510  mode = m_bdata.data();
511  boost::shared_ptr< Points<NekDouble> > m_points = PointsManager()[PointsKey(numModes, eGaussLobattoLegendre)];
512  const Array<OneD, const NekDouble>& zp(m_points->GetZ());
513 
514  for (p=0; p<numModes; ++p, mode += numPoints)
515  {
516  for(q = 0; q < numPoints; ++q)
517  {
518  mode[q] = Polylib::hglj(p, z[q], zp.data(), numModes, 0.0, 0.0);
519  }
520  }
521 
522  // define derivative basis
523  Blas::Dgemm('n', 'n', numPoints, numModes, numPoints, 1.0,
524  D, numPoints, m_bdata.data(), numPoints, 0.0,
525  m_dbdata.data(), numPoints);
526 
527  }//end scope
528  break;
529  case eGauss_Lagrange:
530  {
531  mode = m_bdata.data();
532  boost::shared_ptr< Points<NekDouble> > m_points = PointsManager()[PointsKey(numModes, eGaussGaussLegendre)];
533  const Array<OneD, const NekDouble>& zp(m_points->GetZ());
534 
535  for (p=0; p<numModes; ++p,mode += numPoints)
536  {
537  for(q = 0; q < numPoints; ++q)
538  {
539  mode[q] = Polylib::hgj(p, z[q], zp.data(), numModes, 0.0, 0.0);
540  }
541  }
542 
543  // define derivative basis
544  Blas::Dgemm('n', 'n', numPoints, numModes, numPoints, 1.0,
545  D, numPoints, m_bdata.data(), numPoints, 0.0,
546  m_dbdata.data(), numPoints);
547 
548  }//end scope
549  break;
550  case eFourier:
551 
552  ASSERTL0(numModes%2==0, "Fourier modes should be a factor of 2");
553 
554  for(i = 0; i < numPoints; ++i)
555  {
556  m_bdata[i] = 1.0;
557  m_bdata[numPoints+i] = 0.0;
558 
559  m_dbdata[i] = m_dbdata[numPoints+i] = 0.0;
560  }
561 
562  for (p=1; p < numModes/2; ++p)
563  {
564  for(i = 0; i < numPoints; ++i)
565  {
566  m_bdata[ 2*p *numPoints+i] = cos(p*M_PI*z[i]);
567  m_bdata[(2*p+1)*numPoints+i] = -sin(p*M_PI*z[i]);
568 
569  m_dbdata[ 2*p *numPoints+i] = -p*M_PI*sin(p*M_PI*z[i]);
570  m_dbdata[(2*p+1)*numPoints+i] = -p*M_PI*cos(p*M_PI*z[i]);
571  }
572  }
573 
574  break;
575 
576 
577  // Fourier Single Mode (1st mode)
578  case eFourierSingleMode:
579 
580  for(i = 0; i < numPoints; ++i)
581  {
582  m_bdata[i] = cos(M_PI*z[i]);
583  m_bdata[numPoints+i] = -sin(M_PI*z[i]);
584 
585  m_dbdata[i] = -M_PI*sin(M_PI*z[i]);
586  m_dbdata[numPoints+i] = -M_PI*cos(M_PI*z[i]);
587  }
588 
589  for (p=1; p < numModes/2; ++p)
590  {
591  for(i = 0; i < numPoints; ++i)
592  {
593  m_bdata[ 2*p *numPoints+i] = 0.;
594  m_bdata[(2*p+1)*numPoints+i] = 0.;
595 
596  m_dbdata[ 2*p *numPoints+i] = 0.;
597  m_dbdata[(2*p+1)*numPoints+i] = 0.;
598  }
599  }
600  break;
601 
602  //Fourier Real Half Mode
603  case eFourierHalfModeRe:
604  m_bdata[0] = cos(M_PI*z[0]);
605  m_dbdata[0] = -M_PI*sin(M_PI*z[0]);
606  break;
607 
608  //Fourier Imaginary Half Mode
609  case eFourierHalfModeIm:
610  m_bdata[0] = -sin(M_PI*z[0]);
611  m_dbdata[0] = -M_PI*cos(M_PI*z[0]);
612  break;
613 
614  case eChebyshev:
615  {
616  mode = m_bdata.data();
617 
618  for (p=0,scal = 1; p<numModes; ++p,mode += numPoints)
619  {
620  Polylib::jacobfd(numPoints, z.data(), mode, NULL, p, -0.5, -0.5);
621 
622  for(i = 0; i < numPoints; ++i)
623  {
624  mode[i] *= scal;
625  }
626 
627  scal *= 4*(p+1)*(p+1)/(2*p+2)/(2*p+1);
628  }
629 
630  // Define derivative basis
631  Blas::Dgemm('n', 'n', numPoints, numModes, numPoints, 1.0,
632  D, numPoints, m_bdata.data(), numPoints, 0.0,
633  m_dbdata.data(), numPoints);
634  }
635  break;
636 
637  case eMonomial:
638  {
639  int P = numModes - 1;
640  NekDouble *mode = m_bdata.data();
641 
642  for( int p = 0; p <= P; ++p, mode += numPoints )
643  {
644  for( int i = 0; i < numPoints; ++i )
645  {
646  mode[i] = pow(z[i], p);
647  }
648  }
649 
650  // define derivative basis
651  Blas::Dgemm('n', 'n', numPoints, numModes, numPoints, 1.0,
652  D, numPoints, m_bdata.data(), numPoints, 0.0,
653  m_dbdata.data(),numPoints);
654  }//end scope
655  break;
656  default:
657  ASSERTL0(false, "Basis Type not known or "
658  "not implemented at this time.");
659  }
660  }
661 
662  /** \brief Determine if polynomial basis can be eactly integrated
663  * with itself
664  */
665  bool BasisKey::ExactIprodInt(void) const
666  {
667  bool returnval = false;
668 
669  switch(GetPointsType())
670  {
671  case eGaussGaussLegendre:
672  case eGaussRadauMLegendre:
673  case eGaussRadauPLegendre:
674  case eGaussLobattoLegendre:
675  case eGaussKronrodLegendre:
676  case eGaussRadauKronrodMLegendre:
677  case eGaussLobattoKronrodLegendre:
678  returnval = (GetNumPoints() >= GetNumModes());
679  break;
680 
681  default:
682  break;
683  }
684 
685  return returnval;
686  }
687 
688  /** \brief Determine if basis has collocation property,
689  * i.e. GLL_Lagrange with Lobatto integration of appropriate order,
690  * Gauss_Lagrange with Gauss integration of appropriate order.
691  */
692  bool BasisKey::Collocation() const
693  {
694  return ((m_basistype == eGLL_Lagrange &&
695  GetPointsType() == eGaussLobattoLegendre &&
696  GetNumModes() == GetNumPoints()) ||
697  (m_basistype == eGauss_Lagrange &&
698  GetPointsType() == eGaussGaussLegendre &&
699  GetNumModes() == GetNumPoints()));
700  }
701 
702  // BasisKey compared to BasisKey
703  bool operator == (const BasisKey& x, const BasisKey& y)
704  {
705  return (x.GetPointsKey() == y.GetPointsKey() &&
706  x.m_basistype == y.m_basistype &&
707  x.GetNumModes() == y.GetNumModes());
708  }
709 
710  // BasisKey* compared to BasisKey
711  bool operator == (const BasisKey* x, const BasisKey& y)
712  {
713  return (*x == y);
714  }
715 
716  // \brief BasisKey compared to BasisKey*
717  bool operator == (const BasisKey& x, const BasisKey *y)
718  {
719  return (x == *y);
720  }
721 
722  // \brief BasisKey compared to BasisKey
723  bool operator != (const BasisKey& x, const BasisKey& y)
724  {
725  return (!(x == y));
726  }
727 
728  // BasisKey* compared to BasisKey
729  bool operator != (const BasisKey* x, const BasisKey& y)
730  {
731  return (!(*x == y));
732  }
733 
734  // BasisKey compared to BasisKey*
735  bool operator != (const BasisKey& x, const BasisKey* y)
736  {
737  return (!(x == *y));
738  }
739 
740  } // end of namespace stdregion
741 } // end of namespace stdregion
742 
Generated on Sun Mar 15 2015 00:11:29 for Nektar++ by   doxygen 1.8.1.2