Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
NodalUtil.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File NodalUtil.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: 2D Nodal Triangle Fekete Utilities --
33 // Basis function, Interpolation, Integral, Derivation, etc.
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
37 #include <iomanip>
38 #include <limits>
39 
40 #include <LibUtilities/Foundations/NodalUtil.h>
41 
42 #include <LibUtilities/LinearAlgebra/NekLinSys.hpp>
43 #include <LibUtilities/BasicUtils/SharedArray.hpp>
44 #include <LibUtilities/LinearAlgebra/NekMatrix.hpp>
45 #include <LibUtilities/LinearAlgebra/NekVector.hpp>
46 
47 namespace Nektar
48 {
49  namespace LibUtilities
50  {
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 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 Array<OneD, const 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 Array<OneD, const 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  NekVector<NekDouble> MakeDubinerQuadratureSystem(int nBasisFunctions)
237  {
238  // Make the vector of integrals: note that each modal basis function integrates to zero except for the 0th degree
239  NekVector<NekDouble> g(nBasisFunctions, 0.0);
240  g(0) = sqrt(2.0);
241 
242  return g;
243  }
244 
245  NekVector<NekDouble> MakeTetQuadratureSystem(int nBasisFunctions)
246  {
247  NekVector<NekDouble> g(nBasisFunctions, 0.0);
248  g(0) = 1.0;
249  return g;
250  }
251 
252 
253  NekVector<NekDouble> JacobiPoly(int degree, const NekVector<NekDouble>& x, NekDouble alpha, NekDouble beta)
254  {
255  int size = GetSize(x);
256  NekVector<NekDouble> y(size);
257 
258  if(degree == 0)
259  {
260  // Set y to ones
261  y = NekVector<NekDouble>(size, 1.0);
262 
263  }
264  else if (degree == 1)
265  {
266  for(int el=0; el<size; ++el)
267  {
268  y[el] = 0.5*(alpha - beta + (alpha + beta + 2.0) * x[el]);
269  }
270 
271  }
272  else if (degree > 1)
273  {
274  NekDouble degm1 = degree - 1.0;
275  NekDouble tmp = 2.0 * degm1 + alpha + beta;
276  NekDouble a1 = 2.0 * (degm1 + 1.0) * (degm1 + alpha + beta + 1.0) * tmp;
277  NekDouble a2 = (tmp + 1.0) * (alpha * alpha - beta * beta);
278  NekDouble a3 = tmp * (tmp + 1.0) * (tmp + 2.0);
279  NekDouble a4 = 2.0 * (degm1 + alpha) * (degm1 + beta) * (tmp + 2.0);
280 
281  // TODO: Make this efficient: Reuse the degree-2 for the degree-1 call.
282  // This can be done in linear time, but it is currently implemented to run in exponential time.
283  NekVector<NekDouble> z2 = JacobiPoly(degree-2, x, alpha, beta);
284  NekVector<NekDouble> z1 = JacobiPoly(degree-1, x, alpha, beta);
285 
286  for (int el=0; el<size; ++el)
287  {
288  y[el] = ((a2 + a3 * x[el]) * z1[el] - a4 * z2[el])/a1;
289  }
290 
291  }
292  else
293  {
294  cerr << "Bad degree" << endl;
295  }
296 
297  return y;
298  }
299 
300  NekDouble JacobiPoly(int degree, NekDouble x, NekDouble alpha, NekDouble beta)
301  {
302  NekDouble y = 0.0;
303 
304  if(degree == 0)
305  {
306  y = 1.0;
307  }
308 
309  else if (degree == 1)
310  {
311  y = 0.5*(alpha - beta + (alpha + beta + 2.0) * x);
312  }
313 
314  else if (degree > 1)
315  {
316  NekDouble degm1 = degree - 1.0;
317  NekDouble tmp = 2.0 * degm1 + alpha + beta;
318  NekDouble a1 = 2.0 * (degm1 + 1.0) * (degm1 + alpha + beta + 1.0) * tmp;
319  NekDouble a2 = (tmp + 1.0) * (alpha * alpha - beta * beta);
320  NekDouble a3 = tmp * (tmp + 1.0) * (tmp + 2.0);
321  NekDouble a4 = 2.0 * (degm1 + alpha) * (degm1 + beta) * (tmp + 2.0);
322 
323  // TODO: Make this efficient: Reuse the degree-2 for the degree-1 call.
324  // This can be done in linear time, but it is currently implemented to run in exponential time.
325  NekDouble z2 = JacobiPoly(degree-2, x, alpha, beta);
326  NekDouble z1 = JacobiPoly(degree-1, x, alpha, beta);
327 
328  y = ((a2 + a3 * x) * z1 - a4 * z2)/a1;
329  }
330 
331  else
332  {
333  cerr << "Bad degree" << endl;
334  }
335 
336  return y;
337  }
338 
339 
340  NekVector<NekDouble> LegendrePoly(int degree, const NekVector<NekDouble>& x)
341  {
342  int size = GetSize(x);
343  NekVector<NekDouble> y(size);
344 
345  if(degree > 1)
346  {
347  NekVector<NekDouble> a0(size, 1.0);
348  NekVector<NekDouble> a1 = x;
349  NekVector<NekDouble> a2(size);
350 
351  for(int i=2; i<=degree; ++i)
352  {
353  NekDouble b = NekDouble(2.0*i-1.0)/i;
354  NekDouble c = NekDouble(i-1.0)/i;
355 
356  // multiply each elements in matrix
357  a2 = b * Hadamard(a1, x) - c * a0;
358  a0 = a1;
359  a1 = a2;
360  }
361  y = a2;
362  }
363  else if( degree == 1 )
364  {
365  y = x;
366  }
367  else
368  {
369  y = NekVector<NekDouble>(size, 1.0);
370  }
371  return y;
372  }
373 
374  // Triangle orthonormal basis expansion
375  NekVector<NekDouble> DubinerPoly(int p, int q, const NekVector<NekDouble>& x, const NekVector<NekDouble>& y)
376  {
377  // Get the coordinate transform
378  int size = GetSize(x);
379  NekVector<NekDouble> eta_1(size);
380 
381  // Initialize the horizontal coordinate of the triangle (beta in Barycentric coordinates)
382  for(int el=0; el<size; ++el)
383  {
384  if( y[el] < 1.0 - 1e-16 )
385  {
386  eta_1[el] = 2.0*(1.0 + x[el]) / (1.0 - y[el]) - 1.0;
387  }
388  else
389  {
390  eta_1[el] = -1.0; // When y is close to 1, then we have a removeable singularity
391  }
392  }
393 
394  // Initialize the vertical coordinate of the triangle (gamma in Barycentric coordinates)
395  NekVector<NekDouble> eta_2 = y;
396 
397  // Orthogonal Dubiner polynomial
398  int alpha = 2*p + 1; int beta = 0;
399  NekVector<NekDouble> phi(size);
400  NekVector<NekDouble> psi_a = LegendrePoly(p, eta_1);
401  NekVector<NekDouble> upsilon = JacobiPoly(q, eta_2, alpha, beta);
402  NekDouble w = sqrt((2.0 * p + 1.0) * (p + q + 1.0) / 2.0); // Normalizing Orthonormal weight
403 
404  for(int el=0; el<size; ++el)
405  {
406  NekDouble zeta = pow((1.0 - eta_2[el])/2.0, p);
407  NekDouble psi_b = zeta * upsilon[el];
408  phi[el] = w * psi_a[el] * psi_b;
409  }
410  return phi;
411  }
412 
413 
414  // Tetrahedral orthogonal basis expansion
415  NekVector<NekDouble> TetrahedralBasis(int p, int q, int r, const NekVector<NekDouble>& x,
416  const NekVector<NekDouble>& y, const NekVector<NekDouble>& z)
417  {
418  // Get the coordinate transform
419  int size = GetSize(x);
420  NekVector<NekDouble> eta_1(size), eta_2(size);
421 
422  // Initialize the horizontal coordinate of the Tetrahedral
423  for(int el=0; el<size; ++el)
424  {
425  if( z[el] < -y[el] - numeric_limits<NekDouble>::epsilon() )
426  {
427  eta_1[el] = 2.0*(1.0 + x[el])/(-y[el]-z[el]) - 1.0;
428  }
429  else
430  {
431  eta_1[el] = -1.0;
432  }
433  }
434 
435  // Initialize the coordinate of the Tetrahedral
436  for(int el=0; el<size; ++el)
437  {
438  if( z[el] < 1.0 - numeric_limits<NekDouble>::epsilon())
439  {
440  eta_2[el] = 2.0*(1.0 + y[el]) / (1.0 - z[el]) - 1.0;
441  }
442  else
443  {
444  eta_2[el] = -1.0; // When z is close to 1, then we have a removeable singularity
445  eta_1[el] = -1.0; // When z is close to 1, then we have a removeable singularity
446  }
447  }
448 
449  // Initialize the vertical coordinate of the Tetrahedral
450  NekVector<NekDouble> eta_3 = z;
451 
452  // Orthogonal basis polynomial
453  int alpha = 2*p + 1; int beta = 0;
454  int alpha_r = 2*p + 2*q + 2; int beta_r = 0;
455  NekVector<NekDouble> phi(size);
456  NekVector<NekDouble> psi_a = LegendrePoly(p,eta_1);
457  NekVector<NekDouble> psi_bpq = JacobiPoly(q, eta_2, alpha, beta);
458  NekVector<NekDouble> psi_cpqr = JacobiPoly(r, eta_3, alpha_r, beta_r);
459  NekDouble w = 1.0;
460 
461  for(int el=0; el<size; ++el)
462  {
463  NekDouble zeta_1 = pow((1.0 - eta_2[el])/2.0, p);
464  NekDouble zeta_2 = pow((1.0 - eta_3[el])/2.0, p+q);
465 
466  NekDouble psi_b = zeta_1 * psi_bpq[el];
467  NekDouble psi_c = zeta_2 * psi_cpqr[el];
468 
469  phi[el] = w * psi_a[el] * psi_b * psi_c;
470 
471  }
472 
473  return phi;
474 
475  }
476 
477  NekMatrix<NekDouble> GetTetVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y, const NekVector<NekDouble>& z, int degree)
478  {
479  //cout << "Begin GetTetVandermonde" << endl;
480  int rows = GetSize(x);
481  int cols = (degree + 1) * (degree + 2) * (degree + 3)/6;
482 
483  NekMatrix<NekDouble> matV(rows, cols, 0.0);
484 
485  for(int d=0, n=0; d<=degree; ++d)
486  {
487  for(int p=0; p <= d; ++p)
488  {
489  for(int q=0; q <= d - p; ++q, ++n)
490  {
491  int r = d - p - q;
492  NekVector<NekDouble> columnVector = TetrahedralBasis(p, q, r, x, y, z);
493 
494  // cout << "degree = " << degree << ", (d,p,q,r) = (" << d << ", " << p << ", " << q << ", " << r << ")" << endl;
495  // Set n-th column of V to the TetrahedralBasis vector
496  for(int i=0; i<rows; ++i)
497  {
498  matV(i,n) = columnVector[i];
499  }
500  }
501  }
502  }
503  return matV;
504  }
505 
506 
507  NekMatrix<NekDouble> GetTetVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y, const NekVector<NekDouble>& z)
508  {
509  int rows = GetSize(x);
510  int degree = GetTetDegree(rows);
511  return GetTetVandermonde(x, y, z, degree);
512  }
513 
514 
515  NekMatrix<NekDouble> GetVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y, int degree)
516  {
517  int rows = GetSize(x);
518  int cols = (degree + 1) * (degree + 2) / 2;
519 
520  NekMatrix<NekDouble> matV(rows, cols, 0.0);
521 
522  for(int d=0, n=0; d<=degree; ++d)
523  {
524  for(int p=0; p<=d; ++p, ++n)
525  {
526  int q = d - p;
527  NekVector<NekDouble> columnVector = DubinerPoly(p, q, x, y);
528 
529  // Set n-th column of V to the DubinerPoly vector
530  for(int i=0; i<rows; ++i)
531  {
532  matV(i,n) = columnVector[i];
533  }
534  }
535  }
536 
537  return matV;
538  }
539 
540  NekMatrix<NekDouble> GetVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y)
541  {
542  int rows = GetSize(x);
543  int degree = GetDegree(rows);
544  return GetVandermonde( x, y, degree );
545  }
546 
547  SharedNekMatrixPtr MakeVmatrixOfTet(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y, const NekVector<NekDouble>& z)
548  {
549  // Construct the Vandermonde matrix of Tetrahedron
550  SharedNekMatrixPtr vMat(new NekMatrix<NekDouble>( GetTetVandermonde(x, y, z)));
551  return vMat;
552  }
553 
554  SharedNekMatrixPtr MakeVmatrixOfDubinerPolynomial(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y)
555  {
556  // Construct the Vandermonde matrix of DubinerPolynomials
557  SharedNekMatrixPtr vandermonde(new NekMatrix<NekDouble>( GetVandermonde(x, y) ));
558 
559  return vandermonde;
560  }
561 
562  NekVector<NekDouble> MakeTetWeights(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y, const NekVector<NekDouble>& z)
563  {
564 
565  NekVector<NekDouble> g = MakeTetQuadratureSystem(x.GetRows());
566 
567  SharedNekMatrixPtr matV = MakeVmatrixOfTet(x, y, z);
568 
569  // Initialize the linear system solver
570  LinearSystem matL(matV);
571 
572 
573  // Deduce the quadrature weights
574  NekVector<NekDouble> w = matL.SolveTranspose(g);
575 
576  return w;
577  }
578 
579  NekVector<NekDouble> MakeQuadratureWeights(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y)
580  {
581 
582  NekVector<NekDouble> g = MakeDubinerQuadratureSystem(x.GetRows());
583 
584  SharedNekMatrixPtr matV = MakeVmatrixOfDubinerPolynomial( x, y);
585 
586  // Initialize the linear system solver
587  LinearSystem matL(matV);
588 
589  // Deduce the quadrature weights
590  return matL.SolveTranspose(g);
591 
592  }
593 
594  NekMatrix<NekDouble> GetTetInterpolationMatrix(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y, const NekVector<NekDouble>& z,
595  const NekVector<NekDouble>& xi, const NekVector<NekDouble>& yi, const NekVector<NekDouble>& zi)
596  {
597  int nNodes = GetSize(x);
598  int degree = GetTetDegree(nNodes);
599 
600  NekMatrix<NekDouble> matS = GetTetVandermonde(x, y, z); // Square 'short' matrix
601  NekMatrix<NekDouble> matT = GetTetVandermonde(xi, yi, zi, degree); // Coefficient interpolation matrix (tall)
602 
603  NekMatrix<NekDouble> invertMatrix = Invert(matS);
604 
605  // Get the interpolation matrix
606  return matT*invertMatrix;
607  }
608 
609  NekMatrix<NekDouble> GetInterpolationMatrix(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
610  const NekVector<NekDouble>& xi, const NekVector<NekDouble>& yi)
611  {
612  int nNodes = GetSize(x);
613  int degree = GetDegree(nNodes);
614 
615  NekMatrix<NekDouble> matS = GetVandermonde(x, y); // Square 'short' matrix
616  NekMatrix<NekDouble> matT = GetVandermonde(xi, yi, degree); // Coefficient interpolation matrix (tall)
617 
618  NekMatrix<NekDouble> invertMatrix = Invert(matS);
619 
620  // Get the interpolation matrix
621  return matT*invertMatrix;
622  }
623 
624  NekVector<NekDouble> LegendrePolyDerivative(int degree, const NekVector<NekDouble>& x)
625  {
626  int size = GetSize(x);
627  NekVector<NekDouble> y(size), b3(size), a2(size);
628 
629  if(degree >= 3)
630  {
631  NekVector<NekDouble> b0(size, 0.0), b1(size, 1.0), b2 = 3.0 * x;
632  NekVector<NekDouble> a0(size, 1.0), a1 = x;
633 
634  for(int n=3; n<=degree; ++n)
635  {
636  a2 = ((2.0*n - 3.0)/(n - 1.0)) * Hadamard(x, a1) - (n - 2.0)/(n - 1.0) * a0;
637  a0 = a1;
638  a1 = a2;
639 
640  b3 = (2.0*n - 1.0)/n * (Hadamard(b2, x) + a2) - (n - 1.0)/n * b1;
641  b1 = b2;
642  b2 = b3;
643  }
644  y = b3;
645 
646  }
647  else if(degree == 2)
648  {
649  y = 3.0 * x;
650  }
651  else if(degree == 1)
652  {
653  y = NekVector<NekDouble>(size, 1.0);
654 
655  }
656  else
657  {
658  y = NekVector<NekDouble>(size, 0.0);
659  }
660  return y;
661  }
662 
663 
664  NekVector<NekDouble> DubinerPolyXDerivative(int p, int q, const NekVector<NekDouble>& x, const NekVector<NekDouble>& y)
665  {
666  // Get the coordinate transform
667  int size = GetSize(y);
668  NekVector<NekDouble> eta_1(size);
669  NekVector<NekDouble> psi_x(x.GetRows());
670  if(p>0)
671  {
672  // Initialize the horizontal coordinate of the triangle (beta in Barycentric coordinates)
673  for(int el=0; el<size; ++el)
674  {
675  if(y[el] < 1.0 - 1e-16)
676  {
677  eta_1[el] = 2.0*(1.0 + x[el]) / (1.0 - y[el]) - 1.0;
678  }
679  else
680  {
681  eta_1[el] = -1.0; // When y is close to 1, then we have a removeable singularity
682  }
683  }
684  // Initialize the vertical coordinate of the triangle (gamma in Barycentric coordinates)
685  NekVector<NekDouble> eta_2 = y;
686 
687  // Orthogonal Dubiner polynomial x-derivative
688  int alpha = 2*p + 1; int beta = 0;
689  NekVector<NekDouble> psi_a_x = LegendrePolyDerivative(p, eta_1);
690  NekVector<NekDouble> upsilon = JacobiPoly(q, eta_2, alpha, beta);
691  NekDouble w = sqrt((2.0*p + 1.0) * (p + q + 1.0) / 2.0); // Normalizing Orthonormal weight
692 
693  for(int i=0; i<size; ++i)
694  {
695  NekDouble zeta = pow((1.0 - eta_2[i])/2.0, (p-1.0));
696  NekDouble psi_b = zeta*upsilon[i];
697  psi_x[i] = w * psi_a_x[i] * psi_b;
698  }
699  }
700  else
701  {
702  for(int i=0; i<int(x.GetRows()); ++i)
703  {
704  psi_x[i] = 0.0;
705  }
706  }
707 
708  return psi_x;
709  }
710 
711  NekVector<NekDouble> TetXDerivative(int p, int q, int r, const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
712  const NekVector<NekDouble>& z)
713  {
714  // Get the coordinate transform
715  int size = GetSize(y);
716  NekVector<NekDouble> eta_1(size), eta_2(size);
717  NekVector<NekDouble> psi_x(x.GetRows());
718  if(p > 0){
719  // Initialize the horizontal coordinate of the Tetrahedral (beta in Barycentric coordinate)
720  for(int el=0; el<size; ++el)
721  {
722  if( y[el] < -z[el] - numeric_limits<NekDouble>::epsilon())
723  {
724  eta_1[el] = 2.0*(1.0 + x[el])/(-y[el]-z[el]) - 1.0;
725  } else
726  {
727  eta_1[el] = -1.0;
728  }
729  }
730 
731  // Initialize the coordinate of the Tetrahedral (gamma in Barycentric coordinate)
732  for(int el=0; el<size; ++el)
733  {
734  if( z[el] < 1.0 - numeric_limits<NekDouble>::epsilon())
735  {
736  eta_2[el] = 2.0*(1.0 + y[el]) / (1.0 - z[el]) - 1.0;
737  }
738  else
739  {
740  eta_2[el] = -1.0; // When z is close to 1, then we have a removeable singularity
741  }
742  }
743  // Initialize the vertical coordinate of the Tetrahedral (delta in Barycentric coordinate)
744  NekVector<NekDouble> eta_3 = z;
745 
746  // Orthogonal basis polynomial x-derivative
747  int alpha = 2*p + 1; int beta = 0;
748  int alpha_r = 2*p + 2*q + 2; int beta_r = 0;
749  NekVector<NekDouble> psi_a_x = LegendrePolyDerivative(p,eta_1);
750  NekVector<NekDouble> psi_bpq = JacobiPoly(q, eta_2, alpha, beta);
751  NekVector<NekDouble> psi_cpqr = JacobiPoly(r, eta_3, alpha_r, beta_r);
752 
753  for(int el=0; el<size; ++el)
754  {
755  NekDouble jacobi_b = pow((1.0-eta_2[el])/2.0, p-1.0);
756  NekDouble jacobi_c = pow((1.0-eta_3[el])/2.0,p+q-1.0);
757  NekDouble psi_b = jacobi_b * psi_bpq[el];
758  NekDouble psi_c = jacobi_c * psi_cpqr[el];
759  psi_x[el] = psi_a_x[el] * psi_b * psi_c;
760  }
761  }
762  else
763  {
764  for(int el=0; el<int(x.GetRows()); ++el)
765  {
766  psi_x[el] = 0.0;
767  }
768  }
769  return psi_x;
770  }
771 
772  NekMatrix<NekDouble> GetVandermondeForTetXDerivative(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
773  const NekVector<NekDouble>& z, int degree)
774  {
775  int rows = GetSize(x);
776  int cols = (degree + 1) * (degree + 2) * (degree + 3)/6;
777  NekMatrix<NekDouble> matVx(rows, cols, 0.0);
778 
779  for(int d=0, n=0; d<=degree; ++d)
780  {
781  for(int p=0; p <= d; ++p)
782  {
783  for(int q=0; q <= d - p; ++q, ++n)
784  {
785  int r = d - p - q;
786  NekVector<NekDouble> columnVector = TetXDerivative(p, q, r, x, y, z);
787 
788  // Set n-th column of V to the TetrahedralBasis vector
789  for(int i=0; i<rows; ++i)
790  {
791  matVx(i,n) = columnVector[i];
792  }
793  }
794  }
795  }
796  return matVx;
797  }
798 
799  NekMatrix<NekDouble> GetVandermondeForTetXDerivative(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y, const NekVector<NekDouble>& z)
800  {
801  int rows = GetSize(x);
802  int degree = GetTetDegree(rows);
803  return GetVandermondeForTetXDerivative(x, y, z, degree);
804  }
805 
806  Points<NekDouble>::MatrixSharedPtrType GetTetXDerivativeMatrix(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
807  const NekVector<NekDouble>& z, const NekVector<NekDouble>& xi, const NekVector<NekDouble>& yi, const NekVector<NekDouble>& zi)
808  {
809  int nNodes = GetSize(x);
810  int degree = GetTetDegree(nNodes);
811  NekMatrix<NekDouble> matS = GetTetVandermonde(x, y, z); // Square 'short' matrix
812  NekMatrix<NekDouble> matTx = GetVandermondeForTetXDerivative(xi, yi, zi, degree); // Tall matrix
813 
814  NekMatrix<NekDouble> invertMatrix = Invert(matS);
815 
816  // Get the Derivation matrix
817  return Points<NekDouble>::MatrixSharedPtrType( new NekMatrix<NekDouble>(matTx*invertMatrix) );
818  }
819 
820 
821  NekMatrix<NekDouble> GetVandermondeForXDerivative(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y, int degree)
822  {
823  int rows = GetSize(x);
824  int cols = (degree + 1) * (degree + 2)/2;
825  NekMatrix<NekDouble> matVx(rows, cols, 0.0);
826  for(int d=0, n=0; d<=degree; ++d)
827  {
828  for(int p=0; p<=d; ++p, ++n)
829  {
830  int q = d - p;
831  NekVector<NekDouble> columnVector = DubinerPolyXDerivative(p, q, x, y);
832 
833  // Set n-th column of Vx to the DubinerPolyXDerivative vector
834  for(int i=0; i<rows; ++i)
835  {
836  matVx(i, n) = columnVector[i];
837  }
838  }
839  }
840 
841  return matVx;
842  }
843 
844 
845  NekMatrix<NekDouble> GetVandermondeForXDerivative(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y)
846  {
847  int rows = GetSize(x);
848  int degree = GetDegree(rows);
849  return GetVandermondeForXDerivative(x, y, degree);
850  }
851 
852  Points<NekDouble>::MatrixSharedPtrType GetXDerivativeMatrix(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
853  const NekVector<NekDouble>& xi, const NekVector<NekDouble>& yi)
854  {
855  int nNodes = GetSize(x);
856  int degree = GetDegree(nNodes);
857  NekMatrix<NekDouble> matS = GetVandermonde(x, y); // Square 'short' matrix
858  NekMatrix<NekDouble> matTx = GetVandermondeForXDerivative(xi, yi, degree); // Tall matrix
859 
860  NekMatrix<NekDouble> invertMatrix = Invert(matS);
861 
862  // Get the Derivation matrix
863  return Points<NekDouble>::MatrixSharedPtrType( new NekMatrix<NekDouble>(matTx*invertMatrix) );
864  }
865 
866 
867  NekVector<NekDouble> JacobiPolyDerivative(int degree, const NekVector<NekDouble>& x, int alpha, int beta)
868  {
869  int size = GetSize(x);
870  NekVector<NekDouble> y(size);
871 
872  if(degree == 0)
873  {
874  y = NekVector<NekDouble>(size, 0.0);
875 
876  }
877  else
878  {
879  y = 0.5 * (alpha + beta + degree + 1) * JacobiPoly(degree - 1, x, alpha + 1, beta + 1);
880  }
881  return y;
882  }
883 
884 
885 
886  NekVector<NekDouble> TetYDerivative(int p, int q, int r, const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
887  const NekVector<NekDouble>& z)
888  {
889  // Get the coordinate transform
890  int size = GetSize(y);
891  NekVector<NekDouble> eta_1(size), eta_2(size);
892  NekVector<NekDouble> eta_1_dy(size), eta_2_dy(size);
893 
894  // Initialize the collapsed horizontal coordinate of the Tetrahedral
895  for(int el=0; el<size; ++el)
896  {
897  if( y[el] < -z[el] - numeric_limits<NekDouble>::epsilon())
898  {
899  eta_1[el] = 2.0*(1.0 + x[el])/(-y[el]-z[el]) - 1.0;
900  eta_1_dy[el] = 2.0*(1.0 + x[el])/((y[el]+z[el])*(y[el]+z[el]));
901  }
902  else
903  {
904  eta_1[el] = -1.0;
905  eta_1_dy[el] = 0.0; // No change in the squeeze direction
906  }
907  }
908 
909  // Initialize the collapsed depth coordinate of the Tetrahedral
910  for(int el=0; el<size; ++el)
911  {
912  if( z[el] < 1.0 - numeric_limits<NekDouble>::epsilon())
913  {
914  eta_2[el] = 2.0*(1.0 + y[el]) / (1.0 - z[el]) - 1.0;
915  eta_2_dy[el] = 2.0/(1.0 - z[el]);
916  }
917  else
918  {
919  // When z is close to 1, then we have a removeable singularity
920  eta_2[el] = -1.0;
921  eta_2_dy[el] = 0.0; // No change in the squeeze direction
922  }
923  }
924  // Initialize the collapsed vertical coordinate of the Tetrahedral
925  NekVector<NekDouble> eta_3 = z;
926 
927  // Orthogonal basis expansion polynomials and their y-derivatives for the tetrahedron
928  // phi(vec xi) = psi_a(eta_1(xi)) * psi_b(eta_2(xi)) * psi_c(eta_3(xi))
929  int alpha_b = 2*p + 1; int beta_b = 0;
930  int alpha_c = 2*p + 2*q + 2; int beta_c = 0;
931  NekVector<NekDouble> ones( size, 1.0 );
932 
933  NekVector<NekDouble> jacobi_b = VectorPower( (ones - eta_2)/2.0, p );
934  NekVector<NekDouble> jacobi_c = VectorPower( (ones - eta_3)/2.0, p+q );
935 
936  NekVector<NekDouble> J_q = JacobiPoly(q, eta_2, alpha_b, beta_b);
937  NekVector<NekDouble> J_r = JacobiPoly(r, eta_3, alpha_c, beta_c);
938 
939  NekVector<NekDouble> psi_a = LegendrePoly(p, eta_1);
940  NekVector<NekDouble> LD = LegendrePolyDerivative(p, eta_1);
941  NekVector<NekDouble> JD = JacobiPolyDerivative(q, eta_2, alpha_b, beta_b);
942  NekVector<NekDouble> psi_b = Hadamard( J_q, jacobi_b );
943  NekVector<NekDouble> psi_c = Hadamard( J_r, jacobi_c );
944 
945  // Compute the partials wrt y (psi_c_dy = 0)
946  NekVector<NekDouble> secondComponentOfPsi_b_dy( size, 0.0 );
947 
948  if(p > 0)
949  {
950  for(int i=0; i<size; ++i)
951  {
952  NekDouble jacobi_b_dy = -p / 2.0 * pow( (1.0 - eta_2[i])/2.0, p-1.0 );
953  secondComponentOfPsi_b_dy[i] = J_q[i] * jacobi_b_dy;
954  }
955  }
956  else
957  {
958  for(int i=0; i<size; ++i)
959  {
960  secondComponentOfPsi_b_dy[i] = 0.0;
961  }
962  }
963 
964  NekVector<NekDouble> psi_a_dy = Hadamard( LegendrePolyDerivative(p, eta_1), eta_1_dy );
965  NekVector<NekDouble> psi_b_dy = Hadamard( Hadamard( JacobiPolyDerivative(q, eta_2, alpha_b, beta_b), jacobi_b)
966  + secondComponentOfPsi_b_dy, eta_2_dy );
967 
968  NekVector<NekDouble> psi_dy(size);
969  for(int k=0; k<size; ++k)
970  {
971  psi_dy[k] = (psi_a_dy[k]*psi_b[k]*psi_c[k]) + (psi_a[k]*psi_b_dy[k]*psi_c[k]);
972  }
973 
974  // Fix singularity at z=1
975  for(int k=0; k<size; ++k)
976  {
977  if( z[k] >= 1.0 - numeric_limits<NekDouble>::epsilon() )
978  {
979  if( p + q > 0 )
980  {
981  psi_dy[k] = ((2.0*p+q+2.0)*JacobiPoly(q, eta_2[k], alpha_b, 1) - (p+q+2.0)*J_q[k]) *
982  psi_a[k] * pow( (1.0-eta_3[k])/2.0, p+q-1 ) / 2.0 * ((p+q+r+3.0)*J_r[k] - (2.0*p+2.0*q+r+3.0)*JacobiPoly(r, eta_3[k], alpha_c, 1));
983  }
984  else
985  {
986  psi_dy[k] = 0.0;
987  }
988  }
989  }
990 
991 // cout << "(p,q,r) = (" << p << ", " << q << ", " << r << ")" << endl;
992 //
993 // cout << "psi_a = " << VectorToString(psi_a) << endl;
994 // cout << "psi_b = " << VectorToString(psi_b) << endl;
995 // cout << "psi_c = " << VectorToString(psi_c) << endl;
996 // cout << "psi_a_dy = " << VectorToString(psi_a_dy) << endl;
997 // cout << "psi_b_dy = " << VectorToString(psi_b_dy) << endl;
998 // cout << "(psi_a_dy*psi_b*psi_c) = " << VectorToString(Hadamard(Hadamard(psi_a_dy,psi_b), psi_c)) << endl;
999 // cout << "(psi_a*jacobi_b*psi_b_dy*psi_c) = " << VectorToString(Hadamard(Hadamard(psi_a,jacobi_b),Hadamard(psi_b_dy,psi_c))) << endl;
1000 // cout << "secondComponentOfPsi_b_dy = " << VectorToString(secondComponentOfPsi_b_dy) << endl;
1001 //
1002 // cout << "psi_dy = " << VectorToString(psi_dy) << endl;
1003 // // // cout << "jacobi_c = " << VectorToString(jacobi_c) << endl;
1004 // // // cout << "eta_3 = " << VectorToString(eta_3) << endl;
1005 // // // cout << "ones - eta_3 = " << VectorToString(ones - eta_3) << endl;
1006 // // // cout << "(ones - eta_3)/2.0 = " << VectorToString((ones - eta_3)/2.0) << endl;
1007 
1008 
1009  return psi_dy;
1010  }
1011 
1012 
1013  NekMatrix<NekDouble> GetVandermondeForTetYDerivative(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
1014  const NekVector<NekDouble>& z, int degree){
1015  int rows = GetSize(y);
1016  int cols = (degree + 1) * (degree + 2) * (degree + 3)/6;
1017  NekMatrix<NekDouble> matVy(rows, cols, 0.0);
1018 
1019  for(int d=0, n=0; d<=degree; ++d)
1020  {
1021  for(int p=0; p <= d; ++p)
1022  {
1023  for(int q=0; q <= d - p; ++q, ++n)
1024  {
1025  int r = d - p - q;
1026  NekVector<NekDouble> columnVector = TetYDerivative(p, q, r, x, y, z);
1027 
1028  // Set n-th column of V to the TetrahedralBasis vector
1029  for(int i=0; i<rows; ++i)
1030  {
1031  matVy(i,n) = columnVector[i];
1032  }
1033  }
1034  }
1035  }
1036  return matVy;
1037  }
1038 
1039  NekMatrix<NekDouble> GetVandermondeForTetYDerivative(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
1040  const NekVector<NekDouble>& z)
1041  {
1042  int rows = GetSize(y);
1043  int degree = GetTetDegree(rows);
1044  return GetVandermondeForTetYDerivative(x, y, z, degree);
1045  }
1046 
1047  Points<NekDouble>::MatrixSharedPtrType GetTetYDerivativeMatrix(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
1048  const NekVector<NekDouble>& z, const NekVector<NekDouble>& xi, const NekVector<NekDouble>& yi, const NekVector<NekDouble>& zi)
1049  {
1050  // int nNodes = GetSize(y);
1051 
1052  NekMatrix<NekDouble> matS = GetTetVandermonde(x, y, z); // Square 'short' matrix
1053 
1054  NekMatrix<NekDouble> matTy = GetVandermondeForTetYDerivative(xi, yi, zi); // Tall matrix
1055 
1056  NekMatrix<NekDouble> invertMatrix = Invert(matS);
1057 
1058  NekMatrix<NekDouble> makeDerivativeMatrix = matTy*invertMatrix;
1059 
1060  return Points<NekDouble>::MatrixSharedPtrType(new NekMatrix<NekDouble>(makeDerivativeMatrix));
1061 
1062  }
1063 
1064  NekVector<NekDouble> TetZDerivative(int p, int q, int r, const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
1065  const NekVector<NekDouble>& z)
1066  {
1067  int size = GetSize(z);
1068  NekVector<NekDouble> eta_1(size), eta_2(size);
1069  NekVector<NekDouble> eta_1_dz(size), eta_2_dz(size);
1070 
1071  // Initialize the collapsed horizontal coordinate of the Tetrahedral
1072  for(int el=0; el<size; ++el)
1073  {
1074  if( y[el] < -z[el] - numeric_limits<NekDouble>::epsilon())
1075  {
1076  eta_1[el] = 2.0*(1.0 + x[el]) / (-y[el]-z[el]) - 1.0;
1077  eta_1_dz[el] = 2.0*(1.0 + x[el]) / ((y[el]+z[el])*(y[el]+z[el]));
1078  }
1079  else
1080  {
1081  eta_1[el] = -1.0;
1082  eta_1_dz[el] = 0.0; // No change in the squeeze direction
1083  }
1084  }
1085 
1086  // Initialize the collapsed depth coordinate of the Tetrahedral
1087  for(int el=0; el<size; ++el)
1088  {
1089  if( z[el] < 1.0 - numeric_limits<NekDouble>::epsilon())
1090  {
1091  eta_2[el] = 2.0*(1.0 + y[el]) / (1.0 - z[el]) - 1.0;
1092  eta_2_dz[el] = 2.0*(1.0 + y[el]) / ((1.0 - z[el])*(1.0 - z[el]));
1093  }
1094  else
1095  { // When z is close to 1, then we have a removeable singularity
1096  eta_2[el] = -1.0;
1097  eta_2_dz[el] = 0.0; // No change in the squeeze direction
1098  }
1099  }
1100 
1101 
1102  // Initialize the collapsed vertical coordinate of the Tetrahedral
1103  NekVector<NekDouble> eta_3 = z;
1104  NekVector<NekDouble> eta_3_dz(size, 1.0);
1105 
1106 
1107  // Orthogonal basis expansion polynomials and their z-derivatives for the tetrahedron
1108  // phi(vec xi) = psi_a(eta_1(xi)) * psi_b(eta_2(xi)) * psi_c(eta_3(xi))
1109  int alpha_b = 2*p + 1; int beta_b = 0;
1110  int alpha_c = 2*p + 2*q + 2; int beta_c = 0;
1111  NekVector<NekDouble> ones( size, 1.0 );
1112 
1113  NekVector<NekDouble> jacobi_b = VectorPower( (ones - eta_2)/2.0, p );
1114  NekVector<NekDouble> jacobi_c = VectorPower( (ones - eta_3)/2.0, p+q );
1115 
1116  NekVector<NekDouble> J_q = JacobiPoly(q, eta_2, alpha_b, beta_b);
1117  NekVector<NekDouble> J_r = JacobiPoly(r, eta_3, alpha_c, beta_c);
1118 
1119  NekVector<NekDouble> psi_a = LegendrePoly(p, eta_1);
1120  NekVector<NekDouble> psi_b = Hadamard( J_q, jacobi_b );
1121  NekVector<NekDouble> psi_c = Hadamard( J_r, jacobi_c );
1122 
1123 
1124  // Compute the partials wrt y and z (psi_c_dy = 0)
1125  NekVector<NekDouble> secondComponentOfPsi_b_dz( size, 0.0 );
1126  NekVector<NekDouble> secondComponentOfPsi_c_dz( size, 0.0 );
1127  if(p > 0)
1128  {
1129  for(int i=0; i<size; ++i)
1130  {
1131  NekDouble jacobi_b_dz = -p / 2.0 * pow( (1.0 - eta_2[i])/2.0, p-1.0 );
1132  secondComponentOfPsi_b_dz[i] = J_q[i] * jacobi_b_dz;
1133  }
1134  }
1135 
1136  if(p + q > 0)
1137  {
1138  for(int i=0; i<size; ++i)
1139  {
1140  NekDouble jacobi_c_dz = -(p+q)/2.0*pow( (1.0 - eta_3[i])/2.0, p+q-1.0);
1141  secondComponentOfPsi_c_dz[i] = J_r[i] * jacobi_c_dz;
1142  }
1143  }
1144 
1145  NekVector<NekDouble> psi_a_dz = Hadamard( LegendrePolyDerivative(p, eta_1), eta_1_dz );
1146  NekVector<NekDouble> psi_b_dz = Hadamard( Hadamard(JacobiPolyDerivative(q, eta_2, alpha_b, beta_b), jacobi_b)
1147  + secondComponentOfPsi_b_dz, eta_2_dz );
1148  NekVector<NekDouble> psi_c_dz = Hadamard( Hadamard(JacobiPolyDerivative(r, eta_3, alpha_c, beta_c), jacobi_c)
1149  + secondComponentOfPsi_c_dz, eta_3_dz );
1150 
1151  NekVector<NekDouble> psi_dz(size);
1152 
1153  for(int k=0; k<size; ++k)
1154  {
1155  psi_dz[k] = psi_a_dz[k] * psi_b[k] * psi_c[k] + psi_a[k] * psi_b_dz[k] * psi_c[k] + psi_a[k] * psi_b[k] * psi_c_dz[k];
1156  }
1157 
1158  return psi_dz;
1159  }
1160 
1161  NekMatrix<NekDouble> GetVandermondeForTetZDerivative(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
1162  const NekVector<NekDouble>& z, int degree){
1163  int rows = GetSize(z);
1164  int cols = (degree + 1) * (degree + 2) * (degree + 3)/6;
1165  NekMatrix<NekDouble> matVz(rows, cols, 0.0);
1166 
1167  for(int d=0, n=0; d<=degree; ++d)
1168  {
1169  for(int p=0; p <= d; ++p)
1170  {
1171  for(int q=0; q <= d - p; ++q, ++n)
1172  {
1173  int r = d - p - q;
1174  NekVector<NekDouble> columnVector = TetZDerivative(p, q, r, x, y, z);
1175 
1176  // Set n-th column of V to the TetrahedralBasis vector
1177  for(int i=0; i<rows; ++i)
1178  {
1179  matVz(i,n) = columnVector[i];
1180  }
1181  }
1182  }
1183  }
1184  return matVz;
1185  }
1186 
1187  NekMatrix<NekDouble> GetVandermondeForTetZDerivative(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
1188  const NekVector<NekDouble>& z)
1189  {
1190  int rows = GetSize(z);
1191  int degree = GetTetDegree(rows);
1192  return GetVandermondeForTetZDerivative(x, y, z, degree);
1193  }
1194 
1195  Points<NekDouble>::MatrixSharedPtrType GetTetZDerivativeMatrix(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
1196  const NekVector<NekDouble>& z, const NekVector<NekDouble>& xi, const NekVector<NekDouble>& yi, const NekVector<NekDouble>& zi)
1197  {
1198  int nNodes = GetSize(z);
1199  int degree = GetTetDegree(nNodes);
1200 
1201 
1202  NekMatrix<NekDouble> matS = GetTetVandermonde(x, y, z); // Square 'short' matrix
1203 
1204  NekMatrix<NekDouble> matTz = GetVandermondeForTetZDerivative(xi, yi, zi, degree); // Tall matrix
1205 
1206  NekMatrix<NekDouble> invertMatrix = Invert(matS);
1207 
1208  NekMatrix<NekDouble> makeDerivativeMatrix = matTz*invertMatrix;
1209 
1210  Points<NekDouble>::MatrixSharedPtrType TetZDerivative;
1211  TetZDerivative = Points<NekDouble>::MatrixSharedPtrType(new NekMatrix<NekDouble>(makeDerivativeMatrix));
1212 
1213  return TetZDerivative;
1214 
1215  }
1216 
1217  NekVector<NekDouble> DubinerPolyYDerivative(int p, int q, const NekVector<NekDouble>& x, const NekVector<NekDouble>& y)
1218  {
1219  // Get the coordinate transform
1220  int size = GetSize(y);
1221  NekVector<NekDouble> eta_1(size);
1222 
1223  // Initialize the horizontal coordinate of the triangle (beta in Barycentric coordinates)
1224  for(int el=0; el<size; ++el)
1225  {
1226  if(y[el] < 1.0 - 1e-16)
1227  {
1228  eta_1[el] = 2.0*(1.0 + x[el]) / (1.0 - y[el]) - 1.0;
1229  }
1230  else
1231  {
1232  eta_1[el] = -1.0; // When y is close to 1, then we have a removeable singularity
1233  }
1234  }
1235  // Initialize the vertical coordinate of the triangle (gamma in Barycentric coordinates)
1236  NekVector<NekDouble> eta_2 = y;
1237 
1238  // Orthogonal Dubiner y-polynomial
1239  int alpha = 2*p + 1; int beta = 0;
1240 
1241  NekVector<NekDouble> psi_a = LegendrePoly(p, eta_1);
1242  NekVector<NekDouble> psi_a_y = LegendrePolyDerivative(p, eta_1);
1243  NekVector<NekDouble> psi_b = JacobiPoly(q, eta_2, alpha, beta);
1244  NekVector<NekDouble> psi_b_y = JacobiPolyDerivative(q, eta_2, alpha, beta);
1245  NekVector<NekDouble> secondComponentOf_psi_b(size);
1246  NekVector<NekDouble> psi_y(size);
1247 
1248  NekVector<NekDouble> first_part_derivative(size);
1249 
1250  if(p > 0)
1251  {
1252  for(int i=0; i<size; ++i)
1253  {
1254 
1255  first_part_derivative[i] = (1.0 + eta_1[i])/2.0 * psi_a_y[i] * pow((1.0 - eta_2[i])/2.0, p-1.0) * psi_b[i];
1256  secondComponentOf_psi_b[i] = -p/2.0 * pow(((1.0 - eta_2[i]) / 2.0), p - 1.0) * psi_b[i];
1257 
1258  }
1259  }
1260  else
1261  {
1262  for(int i=0; i<size; ++i)
1263  {
1264  secondComponentOf_psi_b[i] = 0.0;
1265  first_part_derivative[i] = 0.0;
1266  }
1267  }
1268  for(int k=0; k<size; ++k)
1269  {
1270  psi_y[k] = first_part_derivative[k] + psi_a[k] * ( pow((1.0 - eta_2[k])/2.0, p) * psi_b_y[k] + secondComponentOf_psi_b[k] );
1271 
1272  }
1273  NekDouble w = sqrt((2.0*p + 1.0)*(p + q + 1.0)/2.0); // Normalizing Orthonormal weight
1274 
1275  return w * psi_y;
1276  }
1277 
1278 
1279 
1280  NekMatrix<NekDouble> GetVandermondeForYDerivative(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y, int degree)
1281  {
1282  int rows = GetSize(y);
1283  int cols = (degree + 1) * (degree + 2)/2;
1284  NekMatrix<NekDouble> matVy(rows, cols, 0.0);
1285 
1286  for(int d=0, n=0; d<=degree; ++d)
1287  {
1288  for(int p=0; p<=d; ++p, ++n)
1289  {
1290  int q = d - p;
1291  NekVector<NekDouble> columnVector = DubinerPolyYDerivative(p, q, x, y);
1292 
1293  for(int i=0; i<rows; ++i)
1294  {
1295  matVy(i, n) = columnVector[i];
1296  }
1297  }
1298  }
1299  return matVy;
1300  }
1301 
1302 
1303  Points<NekDouble>::MatrixSharedPtrType GetYDerivativeMatrix(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
1304  const NekVector<NekDouble>& xi, const NekVector<NekDouble>& yi)
1305  {
1306  int nNodes = GetSize(y);
1307  int degree = GetDegree(nNodes);
1308 
1309  NekMatrix<NekDouble> matS = GetVandermonde(x, y); // Square 'short' matrix
1310  NekMatrix<NekDouble> matTy = GetVandermondeForYDerivative(xi, yi, degree); // Tall matrix
1311 
1312  NekMatrix<NekDouble> invertMatrix = Invert(matS);
1313 
1314 
1315  // Get the Derivation matrix
1316  return Points<NekDouble>::MatrixSharedPtrType( new NekMatrix<NekDouble>(matTy*invertMatrix) );
1317  }
1318 
1319 
1320  NekMatrix<NekDouble> GetVandermondeForYDerivative(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y)
1321  {
1322  int rows = GetSize(y);
1323  int degree = GetDegree(rows);
1324  return GetVandermondeForYDerivative(x, y, degree);
1325  }
1326 
1327  NekMatrix<NekDouble> GetMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y, int degree)
1328  {
1329  int rows = GetSize(x);
1330  int cols = (degree + 1)*(degree + 2)/2;
1331  NekMatrix<NekDouble>matV(rows,cols, 0.0);
1332 
1333  for(int d=0, n=0; d <= degree; ++d)
1334  {
1335  for(int p=0; p <= d; ++p, ++n)
1336  {
1337  int q = d - p;
1338 
1339  for(int i=0; i<rows; ++i)
1340  {
1341  matV(i, n) = pow(x[i], p) * pow(y[i], q);
1342  }
1343  }
1344  }
1345  return matV;
1346  }
1347 
1348  NekMatrix<NekDouble> GetMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
1349  const NekVector<NekDouble>& z, int degree)
1350  {
1351  int rows = GetSize(x);
1352  int cols = (degree + 1) * (degree + 2) * (degree + 3)/6;
1353 
1354  NekMatrix<NekDouble> matV(rows, cols, 0.0);
1355 
1356  for(int d=0, n=0; d<=degree; ++d)
1357  {
1358  for(int p=0; p <= d; ++p)
1359  {
1360  for(int q=0; q <= d - p; ++q, ++n)
1361  {
1362  int r = d - p - q;
1363 
1364  // Set n-th column of V to the monomial vector
1365  for(int i=0; i<rows; ++i)
1366  {
1367  matV(i,n) = pow(x[i],p) * pow(y[i],q) * pow(z[i],r);
1368  }
1369  }
1370  }
1371  }
1372  return matV;
1373  }
1374 
1375  NekMatrix<NekDouble> GetMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y, const NekVector<NekDouble>& z)
1376  {
1377  int rows = GetSize(x);
1378  int degree = GetTetDegree(rows);
1379  return GetMonomialVandermonde(x, y, z, degree);
1380  }
1381 
1382  NekMatrix<NekDouble> GetMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y)
1383  {
1384  int rows = GetSize(x);
1385  int degree = GetDegree(rows);
1386  return GetMonomialVandermonde(x, y, degree);
1387  }
1388 
1389  NekMatrix<NekDouble> GetXDerivativeOfMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y, int degree)
1390  {
1391  int rows = GetSize(x);
1392  int cols = (degree + 1) * (degree + 2) / 2;
1393  NekMatrix<NekDouble> matVx(rows, cols, 0.0);
1394 
1395  for(int d=0, n=0; d <= degree; ++d)
1396  {
1397  for(int p=0; p <= d; ++p, ++n)
1398  {
1399  int q = d - p;
1400 
1401  if(p > 0)
1402  {
1403  for(int i=0; i<rows; ++i)
1404  {
1405  matVx(i, n) = p * pow(x[i], p-1.0) * pow(y[i],q);
1406  }
1407  }
1408  else
1409  {
1410  for(int j=0; j<rows; ++j)
1411  {
1412  matVx(j, n) = 0.0;
1413  }
1414  }
1415  }
1416  }
1417  return matVx;
1418  }
1419 
1420  NekMatrix<NekDouble> GetXDerivativeOfMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y)
1421  {
1422  int rows = GetSize(x);
1423  int degree = GetDegree(rows);
1424  return GetXDerivativeOfMonomialVandermonde(x, y, degree);
1425  }
1426 
1427  NekMatrix<NekDouble> GetTetXDerivativeOfMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
1428  const NekVector<NekDouble>& z, int degree)
1429  {
1430  int rows = GetSize(x);
1431  int cols = (degree + 1) * (degree + 2) * (degree + 3) / 6;
1432  NekMatrix<NekDouble> matVx(rows, cols, 0.0);
1433  for(int d=0, n=0; d<=degree; ++d)
1434  {
1435  for(int p=0; p <= d; ++p)
1436  {
1437  for(int q=0; q <= d - p; ++q, ++n)
1438  {
1439  int r = d - p - q;
1440  if(p > 0)
1441  {
1442  // Set n-th column of V to the monomial vector
1443  for(int i=0; i<rows; ++i)
1444  {
1445  matVx(i,n) = p * pow(x[i],p-1.0) * pow(y[i],q) * pow(z[i],r);
1446  }
1447  }
1448  else{
1449  for(int j=0; j<rows; ++j)
1450  {
1451  matVx(j, n) = 0.0;
1452  }
1453  }
1454  }
1455  }
1456  }
1457  return matVx;
1458  }
1459 
1460  NekMatrix<NekDouble> GetTetXDerivativeOfMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
1461  const NekVector<NekDouble>& z)
1462  {
1463  int rows = GetSize(x);
1464  int degree = GetTetDegree(rows);
1465  return GetTetXDerivativeOfMonomialVandermonde(x, y, z, degree);
1466  }
1467 
1468  NekMatrix<NekDouble> GetTetYDerivativeOfMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
1469  const NekVector<NekDouble>& z, int degree)
1470  {
1471  int rows = GetSize(x);
1472  int cols = (degree + 1) * (degree + 2) * (degree + 3) / 6;
1473  NekMatrix<NekDouble> matVy(rows, cols, 0.0);
1474  for(int d=0, n=0; d<=degree; ++d)
1475  {
1476  for(int p=0; p <= d; ++p)
1477  {
1478  for(int q=0; q <= d - p; ++q, ++n)
1479  {
1480  int r = d - p - q;
1481  if(q > 0)
1482  {
1483  // Set n-th column of V to the monomial vector
1484  for(int i=0; i<rows; ++i)
1485  {
1486  matVy(i,n) = q * pow(x[i],p) * pow(y[i],q-1.0) * pow(z[i],r);
1487  }
1488  }
1489  else
1490  {
1491  for(int j=0; j<rows; ++j)
1492  {
1493  matVy(j, n) = 0.0;
1494  }
1495  }
1496  }
1497  }
1498  }
1499  return matVy;
1500  }
1501 
1502  NekMatrix<NekDouble> GetTetYDerivativeOfMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
1503  const NekVector<NekDouble>& z)
1504  {
1505  int rows = GetSize(x);
1506  int degree = GetTetDegree(rows);
1507  return GetTetYDerivativeOfMonomialVandermonde(x, y, z, degree);
1508  }
1509 
1510  NekMatrix<NekDouble> GetTetZDerivativeOfMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
1511  const NekVector<NekDouble>& z, int degree)
1512  {
1513  int rows = GetSize(x);
1514  int cols = (degree + 1) * (degree + 2) * (degree + 3) / 6;
1515  NekMatrix<NekDouble> matVz(rows, cols, 0.0);
1516  for(int d=0, n=0; d<=degree; ++d)
1517  {
1518  for(int p=0; p <= d; ++p)
1519  {
1520  for(int q=0; q <= d - p; ++q, ++n)
1521  {
1522  int r = d - p - q;
1523  if(r > 0)
1524  {
1525  // Set n-th column of V to the monomial vector
1526  for(int i=0; i<rows; ++i)
1527  {
1528  matVz(i,n) = r * pow(x[i],p) * pow(y[i],q) * pow(z[i],r-1.0);
1529  }
1530  }
1531  else{
1532  for(int j=0; j<rows; ++j)
1533  {
1534  matVz(j, n) = 0.0;
1535  }
1536  }
1537  }
1538  }
1539  }
1540  return matVz;
1541  }
1542 
1543  NekMatrix<NekDouble> GetTetZDerivativeOfMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y,
1544  const NekVector<NekDouble>& z)
1545  {
1546  int rows = GetSize(x);
1547  int degree = GetTetDegree(rows);
1548  return GetTetZDerivativeOfMonomialVandermonde(x, y, z, degree);
1549  }
1550 
1551  NekMatrix<NekDouble> GetYDerivativeOfMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y, int degree)
1552  {
1553  int rows = GetSize(x);
1554  int cols = (degree + 1) * (degree + 2) / 2;
1555  NekMatrix<NekDouble> matVy(rows, cols, 0.0);
1556 
1557  for(int d=0, n=0; d <= degree; ++d)
1558  {
1559  for(int p=0; p <= d; ++p, ++n)
1560  {
1561  int q = d - p;
1562  if(q > 0)
1563  {
1564  for(int i=0; i<rows; ++i)
1565  {
1566  matVy(i, n) = q * pow(x[i], p) * pow(y[i], q-1.0);
1567  }
1568  }
1569  else
1570  {
1571  for(int j=0; j<rows; ++j)
1572  {
1573  matVy(j, n) = 0.0;
1574  }
1575  }
1576  }
1577  }
1578  return matVy;
1579  }
1580 
1581  NekMatrix<NekDouble> GetYDerivativeOfMonomialVandermonde(const NekVector<NekDouble>& x, const NekVector<NekDouble>& y)
1582  {
1583  int rows = GetSize(x);
1584  int degree = GetDegree(rows);
1585  return GetYDerivativeOfMonomialVandermonde(x, y, degree);
1586  }
1587 
1588  NekVector<NekDouble> GetIntegralOfMonomialVandermonde(int degree)
1589  {
1590  int cols = (degree + 1) * (degree + 2) / 2;
1591  NekVector<NekDouble>integralMVandermonde(cols, 0.0);
1592 
1593  for(int d=0, n=0; d <= degree; ++d)
1594  {
1595  for(int p=0; p <= d; ++p, ++n)
1596  {
1597  int q = d - p;
1598  int sp = 1 - 2 * (p % 2);
1599  int sq = 1 - 2 * (q % 2);
1600 
1601  integralMVandermonde(n) = NekDouble(sp * (p+1) + sq * (q+1) + sp * sq * (p+q+2)) / ((p+1) * (q+1) * (p+q+2));
1602  }
1603  }
1604  return integralMVandermonde;
1605  }
1606  }
1607 }
1608 
Generated on Sun Mar 15 2015 00:11:33 for Nektar++ by   doxygen 1.8.1.2