Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
StdExpUtil.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File StdExpUtil.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:
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <iostream>
37 #include <iomanip>
38 #include <iosfwd>
39 
40 #include <StdRegions/StdExpUtil.h>
41 
42 #include <LibUtilities/LinearAlgebra/NekMatrixFwd.hpp>
43 #include <LibUtilities/LinearAlgebra/Lapack.hpp>
44 #include <LibUtilities/LinearAlgebra/NekLinSys.hpp>
45 #include <LibUtilities/BasicUtils/SharedArray.hpp>
46 
47 
48 namespace Nektar
49 {
50  namespace StdRegions
51  {
52  template< typename T > NekVector<T> GetColumn(const NekMatrix<T> & matA, int n)
53  {
54  NekVector<T> v(matA.GetRows());
55  for( int i=0; i<matA.GetRows(); ++i )
56  {
57  v[i] = matA(i,n);
58  }
59  return v;
60  }
61 
62  NekMatrix<NekDouble> & SetColumn(NekMatrix<NekDouble> & matA, int n, const NekVector<NekDouble> & x)
63  {
64  for( int i=0; i<int(matA.GetRows()); ++i )
65  {
66  matA(i,n) = x[i];
67  }
68  return matA;
69  }
70 
71  // Standard basis vector in R^M
72  NekVector<NekDouble> GetE(int rows, int n)
73  {
74  NekVector<NekDouble> e(rows, 0.0);
75  e(n) = 1;
76  return e;
77  }
78 
79  NekMatrix<NekDouble> Invert(const NekMatrix<NekDouble> & matA)
80  {
81  int rows = matA.GetRows(), columns = matA.GetColumns();
82  NekMatrix<NekDouble> matX(rows,columns);
83 
84  // The linear system solver
85  LinearSystem<NekMatrix<NekDouble> > matL( SharedNekMatrixPtr(new NekMatrix<NekDouble>(matA)) );
86 
87  // Solve each column for the identity matrix
88  for( int j=0; j<columns; ++j )
89  {
90  SetColumn( matX, j, matL.Solve( GetE(rows,j) ) );
91  }
92 
93  return matX;
94  }
95 
96  NekMatrix<NekDouble> GetTranspose(const NekMatrix<NekDouble> & matA)
97  {
98  int rows = matA.GetRows(), columns = matA.GetColumns();
99  NekMatrix<NekDouble> matX(rows,columns);
100 
101  for( int i=0; i<rows; ++i )
102  {
103  for( int j=0; j<columns; ++j )
104  {
105  matX(j,i) = matA(i,j);
106  }
107  }
108  return matX;
109  }
110 
111  int GetSize(const ConstArray<OneD, NekDouble> & x)
112  {
113  return x.num_elements();
114  }
115 
116  int GetSize(const NekVector<NekDouble> & x)
117  {
118  return x.GetRows();
119  }
120 
121  NekVector<NekDouble> ToVector( const ConstArray<OneD, NekDouble> & x )
122  {
123  return NekVector<NekDouble>( GetSize(x), x.data() );
124  }
125 
126  Array<OneD, NekDouble> ToArray( const NekVector<NekDouble> & x )
127  {
128  return Array<OneD, NekDouble>( GetSize(x), x.GetPtr() );
129  }
130 
131  NekVector<NekDouble> Hadamard( const NekVector<NekDouble> & x, const NekVector<NekDouble> & y )
132  {
133  int size = GetSize(x);
134  NekVector<NekDouble> z(size);
135 
136  for( int i=0; i<size; ++i )
137  {
138  z[i] = x[i] * y[i];
139  }
140  return z;
141  }
142 
143  NekVector<NekDouble> VectorPower( const NekVector<NekDouble> & x, NekDouble p )
144  {
145  int size = GetSize(x);
146  NekVector<NekDouble> z(size);
147 
148  for( int i=0; i<size; ++i )
149  {
150  z[i] = pow( x[i], p );
151  }
152 
153  return z;
154  }
155 
156  // Formatted version of matrix ostream output
157  std::string MatrixToString( const NekMatrix<NekDouble> & A, int precision, NekDouble expSigFigs )
158  {
159  stringstream s;
160  s << setprecision(precision);
161  int M = int(A.GetRows()), N = int(A.GetColumns());
162 
163  for(int i=0; i<M; ++i )
164  {
165  for(int j=0; j<N; ++j)
166  {
167  NekDouble a = MakeRound(expSigFigs * A(i, j)) / expSigFigs;
168  s << setw(7) << right << a;
169  if( j < N-1 )
170  {
171  s << ", ";
172  }
173  }
174  if( i < M-1 )
175  {
176  s << "\n";
177  }
178  }
179  return s.str();
180  }
181 
182  // Formatted version of vector ostream output
183  std::string VectorToString( const NekVector<NekDouble> & v, int precision, NekDouble expSigFigs )
184  {
185  stringstream s;
186  s << setprecision(precision) << "[ ";
187  int N = int(v.GetRows());
188  for(int j=0; j<N; ++j )
189  {
190  NekDouble x = MakeRound(expSigFigs * v(j)) / expSigFigs;
191  s << setw(7) << right << x;
192  if( j < N-1 )
193  {
194  s << ", ";
195  }
196  }
197  s << " ]";
198  return s.str();
199  }
200 
201 
202  int GetTriNumPoints(int degree)
203  {
204  return (degree+1) * (degree+2) / 2;
205  }
206 
207  int GetDegree(int nBasisFunctions)
208  {
209  int degree = int(MakeRound((-3.0 + sqrt(1.0 + 8.0*nBasisFunctions))/2.0));
210 
211  // TODO: Find out why ASSERTL0 and ASSERTL1 don't work
212  ASSERTL1( GetTriNumPoints(degree) == nBasisFunctions, "The number of points defines an expansion of fractional degree, which is not supported." );
213  return degree;
214  }
215 
216  int GetTetNumPoints(int degree)
217  {
218  return (degree+1) * (degree+2) * (degree+3) / 6;
219  }
220 
221  // Get Tetrahedral number, where Tn = (d+1)(d+2)(d+3)/6
222  int GetTetDegree(int nBasisFunc)
223  {
224  NekDouble eq = pow( 81.0 * nBasisFunc + 3.0 * sqrt(-3.0 + 729.0 * nBasisFunc * nBasisFunc), 1.0/3.0);
225  int degree = int(MakeRound(eq/3.0 + 1.0/eq - 1.0)) - 1;
226 
227  ASSERTL1( GetTetNumPoints(degree) == nBasisFunc, "The number of points defines an expansion of fractional degree, which is not supported." );
228  return degree;
229  }
230 
231  NekDouble MakeRound(NekDouble x)
232  {
233  return floor(x + 0.5);
234  }
235 
236  NekMatrix<NekDouble> GetMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y, int degree)
237  {
238  int rows = GetSize(x);
239  int cols = (degree + 1)*(degree + 2)/2;
240  NekMatrix<NekDouble>matV(rows,cols, 0.0);
241 
242  for(int d=0, n=0; d <= degree; ++d)
243  {
244  for(int p=0; p <= d; ++p, ++n)
245  {
246  int q = d - p;
247 
248  for(int i=0; i<rows; ++i)
249  {
250  matV(i, n) = pow(x[i], p) * pow(y[i], q);
251  }
252  }
253  }
254  return matV;
255  }
256 
257  NekMatrix<NekDouble> GetMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
258  const NekVector<NekDouble>& z, int degree)
259  {
260  int rows = GetSize(x);
261  int cols = (degree + 1) * (degree + 2) * (degree + 3)/6;
262 
263  NekMatrix<NekDouble> matV(rows, cols, 0.0);
264 
265  for(int d=0, n=0; d<=degree; ++d)
266  {
267  for(int p=0; p <= d; ++p)
268  {
269  for(int q=0; q <= d - p; ++q, ++n)
270  {
271  int r = d - p - q;
272 
273  // Set n-th column of V to the monomial vector
274  for(int i=0; i<rows; ++i)
275  {
276  matV(i,n) = pow(x[i],p) * pow(y[i],q) * pow(z[i],r);
277  }
278  }
279  }
280  }
281  return matV;
282  }
283 
284  NekMatrix<NekDouble> GetMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y, const NekVector<NekDouble>& z)
285  {
286  int rows = GetSize(x);
287  int degree = GetTetDegree(rows);
288  return GetMonomialVandermonde(x, y, z, degree);
289  }
290 
291  NekMatrix<NekDouble> GetMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y)
292  {
293  int rows = GetSize(x);
294  int degree = GetDegree(rows);
295  return GetMonomialVandermonde(x, y, degree);
296  }
297 
298  NekMatrix<NekDouble> GetXDerivativeOfMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y, int degree)
299  {
300  int rows = GetSize(x);
301  int cols = (degree + 1) * (degree + 2) / 2;
302  NekMatrix<NekDouble> matVx(rows, cols, 0.0);
303 
304  for(int d=0, n=0; d <= degree; ++d)
305  {
306  for(int p=0; p <= d; ++p, ++n)
307  {
308  int q = d - p;
309 
310  if(p > 0)
311  {
312  for(int i=0; i<rows; ++i)
313  {
314  matVx(i, n) = p * pow(x[i], p-1.0) * pow(y[i],q);
315  }
316  }
317  else
318  {
319  for(int j=0; j<rows; ++j)
320  {
321  matVx(j, n) = 0.0;
322  }
323  }
324  }
325  }
326  return matVx;
327  }
328 
329  NekMatrix<NekDouble> GetXDerivativeOfMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y)
330  {
331  int rows = GetSize(x);
332  int degree = GetDegree(rows);
333  return GetXDerivativeOfMonomialVandermonde(x, y, degree);
334  }
335 
336  NekMatrix<NekDouble> GetTetXDerivativeOfMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
337  const NekVector<NekDouble>& z, int degree)
338  {
339  int rows = GetSize(x);
340  int cols = (degree + 1) * (degree + 2) * (degree + 3) / 6;
341  NekMatrix<NekDouble> matVx(rows, cols, 0.0);
342  for(int d=0, n=0; d<=degree; ++d)
343  {
344  for(int p=0; p <= d; ++p)
345  {
346  for(int q=0; q <= d - p; ++q, ++n)
347  {
348  int r = d - p - q;
349  if(p > 0)
350  {
351  // Set n-th column of V to the monomial vector
352  for(int i=0; i<rows; ++i)
353  {
354  matVx(i,n) = p * pow(x[i],p-1.0) * pow(y[i],q) * pow(z[i],r);
355  }
356  }
357  else{
358  for(int j=0; j<rows; ++j)
359  {
360  matVx(j, n) = 0.0;
361  }
362  }
363  }
364  }
365  }
366  return matVx;
367  }
368 
369  NekMatrix<NekDouble> GetTetXDerivativeOfMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
370  const NekVector<NekDouble>& z)
371  {
372  int rows = GetSize(x);
373  int degree = GetTetDegree(rows);
374  return GetTetXDerivativeOfMonomialVandermonde(x, y, z, degree);
375  }
376 
377  NekMatrix<NekDouble> GetTetYDerivativeOfMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
378  const NekVector<NekDouble>& z, int degree)
379  {
380  int rows = GetSize(x);
381  int cols = (degree + 1) * (degree + 2) * (degree + 3) / 6;
382  NekMatrix<NekDouble> matVy(rows, cols, 0.0);
383  for(int d=0, n=0; d<=degree; ++d)
384  {
385  for(int p=0; p <= d; ++p)
386  {
387  for(int q=0; q <= d - p; ++q, ++n)
388  {
389  int r = d - p - q;
390  if(q > 0)
391  {
392  // Set n-th column of V to the monomial vector
393  for(int i=0; i<rows; ++i)
394  {
395  matVy(i,n) = q * pow(x[i],p) * pow(y[i],q-1.0) * pow(z[i],r);
396  }
397  }
398  else
399  {
400  for(int j=0; j<rows; ++j)
401  {
402  matVy(j, n) = 0.0;
403  }
404  }
405  }
406  }
407  }
408  return matVy;
409  }
410 
411  NekMatrix<NekDouble> GetTetYDerivativeOfMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
412  const NekVector<NekDouble>& z)
413  {
414  int rows = GetSize(x);
415  int degree = GetTetDegree(rows);
416  return GetTetYDerivativeOfMonomialVandermonde(x, y, z, degree);
417  }
418 
419  NekMatrix<NekDouble> GetTetZDerivativeOfMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
420  const NekVector<NekDouble>& z, int degree)
421  {
422  int rows = GetSize(x);
423  int cols = (degree + 1) * (degree + 2) * (degree + 3) / 6;
424  NekMatrix<NekDouble> matVz(rows, cols, 0.0);
425  for(int d=0, n=0; d<=degree; ++d)
426  {
427  for(int p=0; p <= d; ++p)
428  {
429  for(int q=0; q <= d - p; ++q, ++n)
430  {
431  int r = d - p - q;
432  if(r > 0)
433  {
434  // Set n-th column of V to the monomial vector
435  for(int i=0; i<rows; ++i)
436  {
437  matVz(i,n) = r * pow(x[i],p) * pow(y[i],q) * pow(z[i],r-1.0);
438  }
439  }
440  else{
441  for(int j=0; j<rows; ++j)
442  {
443  matVz(j, n) = 0.0;
444  }
445  }
446  }
447  }
448  }
449  return matVz;
450  }
451 
452  NekMatrix<NekDouble> GetTetZDerivativeOfMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
453  const NekVector<NekDouble>& z)
454  {
455  int rows = GetSize(x);
456  int degree = GetTetDegree(rows);
457  return GetTetZDerivativeOfMonomialVandermonde(x, y, z, degree);
458  }
459 
460  NekMatrix<NekDouble> GetYDerivativeOfMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y, int degree)
461  {
462  int rows = GetSize(x);
463  int cols = (degree + 1) * (degree + 2) / 2;
464  NekMatrix<NekDouble> matVy(rows, cols, 0.0);
465 
466  for(int d=0, n=0; d <= degree; ++d)
467  {
468  for(int p=0; p <= d; ++p, ++n)
469  {
470  int q = d - p;
471  if(q > 0)
472  {
473  for(int i=0; i<rows; ++i)
474  {
475  matVy(i, n) = q * pow(x[i], p) * pow(y[i], q-1.0);
476  }
477  }
478  else
479  {
480  for(int j=0; j<rows; ++j)
481  {
482  matVy(j, n) = 0.0;
483  }
484  }
485  }
486  }
487  return matVy;
488  }
489 
490  NekMatrix<NekDouble> GetYDerivativeOfMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y)
491  {
492  int rows = GetSize(x);
493  int degree = GetDegree(rows);
494  return GetYDerivativeOfMonomialVandermonde(x, y, degree);
495  }
496 
497  NekVector<NekDouble> GetIntegralOfMonomialVandermonde(int degree)
498  {
499  int cols = (degree + 1) * (degree + 2) / 2;
500  NekVector<NekDouble>integralMVandermonde(cols, 0.0);
501 
502  for(int d=0, n=0; d <= degree; ++d)
503  {
504  for(int p=0; p <= d; ++p, ++n)
505  {
506  int q = d - p;
507  int sp = 1 - 2 * (p % 2);
508  int sq = 1 - 2 * (q % 2);
509 
510  integralMVandermonde(n) = NekDouble(sp * (p+1) + sq * (q+1) + sp * sq * (p+q+2)) / ((p+1) * (q+1) * (p+q+2));
511  }
512  }
513  return integralMVandermonde;
514  }
515 
516 
517 
518  } // end of namespace Stdregions
519 } // end of namespace Nektar
Generated on Sun Mar 15 2015 00:11:34 for Nektar++ by   doxygen 1.8.1.2