Nektar++
Functions
Polylib Namespace Reference

The namespace associated with the the Polylib library (Polylib introduction) More...

Functions

double optdiff (double xl, double xr)
 The following function is used to circumvent/reduce "Subtractive Cancellation" The expression 1/dz is replaced by optinvsub(.,.) Added on 26 April 2017. More...
 
double laginterp (double z, int j, const double *zj, int np)
 
static void Jacobz (const int n, double *z, const double alpha, const double beta)
 Calculate the n zeros, z, of the Jacobi polynomial, i.e. \( P_n^{\alpha,\beta}(z) = 0 \). More...
 
static void TriQL (const int n, double *d, double *e, double **z)
 QL algorithm for symmetric tridiagonal matrix. More...
 
double gammaF (const double x)
 Calculate the Gamma function , \( \Gamma(n)\), for integer values and halves. More...
 
double gammaFracGammaF (const int x, const double alpha, const int y, const double beta)
 Calculate fraction of two Gamma functions, \( \Gamma(x+\alpha)/\Gamma(y+\beta) \), for integer values and halves. More...
 
static void RecCoeff (const int n, double *a, double *b, const double alpha, const double beta)
 The routine finds the recurrence coefficients a and b of the orthogonal polynomials. More...
 
void JKMatrix (int n, double *a, double *b)
 Calcualtes the Jacobi-kronrod matrix by determining the a and coefficients. More...
 
void chri1 (int n, double *a, double *b, double *a0, double *b0, double z)
 Given a weight function \(w(t)\) through the first n+1 coefficients a and b of its orthogonal polynomials this routine generates the first n recurrence coefficients for the orthogonal polynomials relative to the modified weight function \((t-z)w(t)\). More...
 
void zwgj (double *z, double *w, const int np, const double alpha, const double beta)
 Gauss-Jacobi zeros and weights. More...
 
void zwgrjm (double *z, double *w, const int np, const double alpha, const double beta)
 Gauss-Radau-Jacobi zeros and weights with end point at z=-1. More...
 
void zwgrjp (double *z, double *w, const int np, const double alpha, const double beta)
 Gauss-Radau-Jacobi zeros and weights with end point at z=1. More...
 
void zwglj (double *z, double *w, const int np, const double alpha, const double beta)
 Gauss-Lobatto-Jacobi zeros and weights with end point at z=-1,1. More...
 
void zwgk (double *z, double *w, const int npt, const double alpha, const double beta)
 Gauss-Kronrod-Jacobi zeros and weights. More...
 
void zwrk (double *z, double *w, const int npt, const double alpha, const double beta)
 Gauss-Radau-Kronrod-Jacobi zeros and weights. More...
 
void zwlk (double *z, double *w, const int npt, const double alpha, const double beta)
 Gauss-Lobatto-Kronrod-Jacobi zeros and weights. More...
 
void Qg (double *Q, const double *z, const int np, const int offset)
 Compute the Integration Matrix. More...
 
void Dgj (double *D, const double *z, const int np, const double alpha, const double beta)
 Compute the Derivative Matrix and its transpose associated with the Gauss-Jacobi zeros. More...
 
void Dgrjm (double *D, const double *z, const int np, const double alpha, const double beta)
 Compute the Derivative Matrix and its transpose associated with the Gauss-Radau-Jacobi zeros with a zero at z=-1. More...
 
void Dgrjp (double *D, const double *z, const int np, const double alpha, const double beta)
 Compute the Derivative Matrix associated with the Gauss-Radau-Jacobi zeros with a zero at z=1. More...
 
void Dglj (double *D, const double *z, const int np, const double alpha, const double beta)
 Compute the Derivative Matrix associated with the Gauss-Lobatto-Jacobi zeros. More...
 
double hgj (const int i, const double z, const double *zgj, const int np, const double alpha, const double beta)
 Compute the value of the i th Lagrangian interpolant through the np Gauss-Jacobi points zgj at the arbitrary location z. More...
 
double hgrjm (const int i, const double z, const double *zgrj, const int np, const double alpha, const double beta)
 Compute the value of the i th Lagrangian interpolant through the np Gauss-Radau-Jacobi points zgrj at the arbitrary location z. This routine assumes zgrj includes the point -1. More...
 
double hgrjp (const int i, const double z, const double *zgrj, const int np, const double alpha, const double beta)
 Compute the value of the i th Lagrangian interpolant through the np Gauss-Radau-Jacobi points zgrj at the arbitrary location z. This routine assumes zgrj includes the point +1. More...
 
double hglj (const int i, const double z, const double *zglj, const int np, const double alpha, const double beta)
 Compute the value of the i th Lagrangian interpolant through the np Gauss-Lobatto-Jacobi points zgrj at the arbitrary location z. More...
 
void Imgj (double *im, const double *zgj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
 Interpolation Operator from Gauss-Jacobi points to an arbitrary distribution at points zm. More...
 
void Imgrjm (double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
 Interpolation Operator from Gauss-Radau-Jacobi points (including z=-1) to an arbitrary distrubtion at points zm. More...
 
void Imgrjp (double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
 Interpolation Operator from Gauss-Radau-Jacobi points (including z=1) to an arbitrary distrubtion at points zm. More...
 
void Imglj (double *im, const double *zglj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
 Interpolation Operator from Gauss-Lobatto-Jacobi points to an arbitrary distrubtion at points zm. More...
 
void polycoeffs (const int i, const double *z, double *c, const int np)
 
void jacobfd (const int np, const double *z, double *poly_in, double *polyd, const int n, const double alpha, const double beta)
 Routine to calculate Jacobi polynomials, \( P^{\alpha,\beta}_n(z) \), and their first derivative, \( \frac{d}{dz} P^{\alpha,\beta}_n(z) \). More...
 
void jacobd (const int np, const double *z, double *polyd, const int n, const double alpha, const double beta)
 Calculate the derivative of Jacobi polynomials. More...
 
void JacZeros (const int n, double *a, double *b, const double alpha, const double beta)
 Zero and Weight determination through the eigenvalues and eigenvectors of a tridiagonal matrix from the three term recurrence relationship. More...
 
std::complex< Nektar::NekDouble > ImagBesselComp (int n, std::complex< Nektar::NekDouble > y)
 Calcualte the bessel function of the first kind with complex double input y. Taken from Numerical Recipies in C. More...
 

Detailed Description

The namespace associated with the the Polylib library (Polylib introduction)

Function Documentation

◆ chri1()

void Polylib::chri1 ( int  n,
double *  a,
double *  b,
double *  a0,
double *  b0,
double  z 
)

Given a weight function \(w(t)\) through the first n+1 coefficients a and b of its orthogonal polynomials this routine generates the first n recurrence coefficients for the orthogonal polynomials relative to the modified weight function \((t-z)w(t)\).

The result will be placed in the array a0 and b0.

Definition at line 1797 of file Polylib.cpp.

1798 {
1799 
1800  double q = ceil(3.0 * n / 2);
1801  int size = (int)q + 1;
1802  if (size < n + 1)
1803  {
1804  fprintf(stderr, "input arrays a and b are too short\n");
1805  }
1806  double *r = new double[n + 1];
1807  r[0] = z - a[0];
1808  r[1] = z - a[1] - b[1] / r[0];
1809  a0[0] = a[1] + r[1] - r[0];
1810  b0[0] = -r[0] * b[0];
1811 
1812  if (n == 1)
1813  {
1814  delete[] r;
1815  return;
1816  }
1817  int k = 0;
1818  for (k = 1; k < n; k++)
1819  {
1820  r[k + 1] = z - a[k + 1] - b[k + 1] / r[k];
1821  a0[k] = a[k + 1] + r[k + 1] - r[k];
1822  b0[k] = b[k] * r[k] / r[k - 1];
1823  }
1824  delete[] r;
1825 }

Referenced by zwlk(), and zwrk().

◆ Dgj()

void Polylib::Dgj ( double *  D,
const double *  z,
const int  np,
const double  alpha,
const double  beta 
)

Compute the Derivative Matrix and its transpose associated with the Gauss-Jacobi zeros.

  • Compute the derivative matrix, d, and its transpose, dt, associated with the n_th order Lagrangian interpolants through the np Gauss-Jacobi points z such that
    \( \frac{du}{dz}(z[i]) = \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \)

Definition at line 627 of file Polylib.cpp.

629 {
630 
631  double one = 1.0, two = 2.0;
632 
633  if (np <= 0)
634  {
635  D[0] = 0.0;
636  }
637  else
638  {
639  int i, j;
640  double *pd;
641 
642  pd = (double *)malloc(np * sizeof(double));
643  jacobd(np, z, pd, np, alpha, beta);
644 
645  for (i = 0; i < np; i++)
646  {
647  for (j = 0; j < np; j++)
648  {
649 
650  if (i != j)
651  D[i * np + j] = pd[j] / (pd[i] * (z[j] - z[i]));
652  else
653  D[i * np + j] =
654  (alpha - beta + (alpha + beta + two) * z[j]) /
655  (two * (one - z[j] * z[j]));
656  }
657  }
658  free(pd);
659  }
660  return;
661 }
Nektar::LibUtilities::beta
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:61
Polylib::jacobd
void jacobd(const int np, const double *z, double *polyd, const int n, const double alpha, const double beta)
Calculate the derivative of Jacobi polynomials.
Definition: Polylib.cpp:1293

References Nektar::LibUtilities::beta, and jacobd().

Referenced by main(), and Nektar::LibUtilities::GaussPoints::v_CalculateDerivMatrix().

◆ Dglj()

void Polylib::Dglj ( double *  D,
const double *  z,
const int  np,
const double  alpha,
const double  beta 
)

Compute the Derivative Matrix associated with the Gauss-Lobatto-Jacobi zeros.

  • Compute the derivative matrix, d, associated with the n_th order Lagrange interpolants through the np Gauss-Lobatto-Jacobi points z such that
    \( \frac{du}{dz}(z[i]) = \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \)

Definition at line 787 of file Polylib.cpp.

789 {
790  if (np <= 1)
791  {
792  D[0] = 0.0;
793  }
794  else
795  {
796  int i, j;
797  double one = 1.0, two = 2.0;
798  double *pd;
799 
800  pd = (double *)malloc(np * sizeof(double));
801 
802  pd[0] = two * pow(-one, np) * gammaFracGammaF(np, beta, np - 1, 0.0);
803  pd[0] /= gammaF(beta + two);
804  jacobd(np - 2, z + 1, pd + 1, np - 2, alpha + 1, beta + 1);
805  for (i = 1; i < np - 1; ++i)
806  pd[i] *= (one - z[i] * z[i]);
807  pd[np - 1] = -two * gammaFracGammaF(np, alpha, np - 1, 0.0);
808  pd[np - 1] /= gammaF(alpha + two);
809 
810  for (i = 0; i < np; i++)
811  {
812  for (j = 0; j < np; j++)
813  {
814  if (i != j)
815  D[i * np + j] = pd[j] / (pd[i] * (z[j] - z[i]));
816  else
817  {
818  if (j == 0)
819  D[i * np + j] =
820  (alpha - (np - 1) * (np + alpha + beta)) /
821  (two * (beta + two));
822  else if (j == np - 1)
823  D[i * np + j] =
824  -(beta - (np - 1) * (np + alpha + beta)) /
825  (two * (alpha + two));
826  else
827  D[i * np + j] = (alpha - beta + (alpha + beta) * z[j]) /
828  (two * (one - z[j] * z[j]));
829  }
830  }
831  }
832  free(pd);
833  }
834 
835  return;
836 }
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer values and halves.
Definition: Polylib.cpp:1322
Polylib::gammaFracGammaF
double gammaFracGammaF(const int, const double, const int, const double)
Calculate fraction of two Gamma functions, , for integer values and halves.
Definition: Polylib.cpp:1371

References Nektar::LibUtilities::beta, gammaF(), gammaFracGammaF(), and jacobd().

Referenced by main(), and Nektar::LibUtilities::GaussPoints::v_CalculateDerivMatrix().

◆ Dgrjm()

void Polylib::Dgrjm ( double *  D,
const double *  z,
const int  np,
const double  alpha,
const double  beta 
)

Compute the Derivative Matrix and its transpose associated with the Gauss-Radau-Jacobi zeros with a zero at z=-1.

  • Compute the derivative matrix, d, associated with the n_th order Lagrangian interpolants through the np Gauss-Radau-Jacobi points z such that
    \( \frac{du}{dz}(z[i]) = \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \)

Definition at line 674 of file Polylib.cpp.

676 {
677 
678  if (np <= 0)
679  {
680  D[0] = 0.0;
681  }
682  else
683  {
684  int i, j;
685  double one = 1.0, two = 2.0;
686  double *pd;
687 
688  pd = (double *)malloc(np * sizeof(double));
689 
690  pd[0] = pow(-one, np - 1) * gammaFracGammaF(np + 1, beta, np, 0.0);
691  pd[0] /= gammaF(beta + two);
692  jacobd(np - 1, z + 1, pd + 1, np - 1, alpha, beta + 1);
693  for (i = 1; i < np; ++i)
694  pd[i] *= (1 + z[i]);
695 
696  for (i = 0; i < np; i++)
697  {
698  for (j = 0; j < np; j++)
699  {
700  if (i != j)
701  D[i * np + j] = pd[j] / (pd[i] * (z[j] - z[i]));
702  else
703  {
704  if (j == 0)
705  D[i * np + j] = -(np + alpha + beta + one) *
706  (np - one) / (two * (beta + two));
707  else
708  D[i * np + j] =
709  (alpha - beta + one + (alpha + beta + one) * z[j]) /
710  (two * (one - z[j] * z[j]));
711  }
712  }
713  }
714  free(pd);
715  }
716 
717  return;
718 }

References Nektar::LibUtilities::beta, gammaF(), gammaFracGammaF(), and jacobd().

Referenced by main(), and Nektar::LibUtilities::GaussPoints::v_CalculateDerivMatrix().

◆ Dgrjp()

void Polylib::Dgrjp ( double *  D,
const double *  z,
const int  np,
const double  alpha,
const double  beta 
)

Compute the Derivative Matrix associated with the Gauss-Radau-Jacobi zeros with a zero at z=1.

  • Compute the derivative matrix, d, associated with the n_th order Lagrangian interpolants through the np Gauss-Radau-Jacobi points z such that
    \( \frac{du}{dz}(z[i]) = \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \)

Definition at line 730 of file Polylib.cpp.

732 {
733 
734  if (np <= 0)
735  {
736  D[0] = 0.0;
737  }
738  else
739  {
740  int i, j;
741  double one = 1.0, two = 2.0;
742  double *pd;
743 
744  pd = (double *)malloc(np * sizeof(double));
745 
746  jacobd(np - 1, z, pd, np - 1, alpha + 1, beta);
747  for (i = 0; i < np - 1; ++i)
748  pd[i] *= (1 - z[i]);
749  pd[np - 1] = -gammaFracGammaF(np + 1, alpha, np, 0.0);
750  pd[np - 1] /= gammaF(alpha + two);
751 
752  for (i = 0; i < np; i++)
753  {
754  for (j = 0; j < np; j++)
755  {
756  if (i != j)
757  D[i * np + j] = pd[j] / (pd[i] * (z[j] - z[i]));
758  else
759  {
760  if (j == np - 1)
761  D[i * np + j] = (np + alpha + beta + one) * (np - one) /
762  (two * (alpha + two));
763  else
764  D[i * np + j] =
765  (alpha - beta - one + (alpha + beta + one) * z[j]) /
766  (two * (one - z[j] * z[j]));
767  }
768  }
769  }
770  free(pd);
771  }
772 
773  return;
774 }

References Nektar::LibUtilities::beta, gammaF(), gammaFracGammaF(), and jacobd().

Referenced by main(), and Nektar::LibUtilities::GaussPoints::v_CalculateDerivMatrix().

◆ gammaF()

double Polylib::gammaF ( const double  )

Calculate the Gamma function , \( \Gamma(n)\), for integer values and halves.

Determine the value of \(\Gamma(n)\) using:

\( \Gamma(n) = (n-1)! \mbox{ or } \Gamma(n+1/2) = (n-1/2)\Gamma(n-1/2)\)

where \( \Gamma(1/2) = \sqrt(\pi)\)

Definition at line 1322 of file Polylib.cpp.

1323 {
1324  double gamma = 1.0;
1325 
1326  if (x == -0.5)
1327  gamma = -2.0 * sqrt(M_PI);
1328  else if (!x)
1329  return gamma;
1330  else if ((x - (int)x) == 0.5)
1331  {
1332  int n = (int)x;
1333  double tmp = x;
1334 
1335  gamma = sqrt(M_PI);
1336  while (n--)
1337  {
1338  tmp -= 1.0;
1339  gamma *= tmp;
1340  }
1341  }
1342  else if ((x - (int)x) == 0.0)
1343  {
1344  int n = (int)x;
1345  double tmp = x;
1346 
1347  while (--n)
1348  {
1349  tmp -= 1.0;
1350  gamma *= tmp;
1351  }
1352  }
1353  else
1354  fprintf(stderr, "%lf is not of integer or half order\n", x);
1355  return gamma;
1356 }
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294

References tinysimd::sqrt().

Referenced by Dglj(), Dgrjm(), Dgrjp(), RecCoeff(), and test_gamma_fraction().

◆ gammaFracGammaF()

double Polylib::gammaFracGammaF ( const int  ,
const double  ,
const int  ,
const double   
)

Calculate fraction of two Gamma functions, \( \Gamma(x+\alpha)/\Gamma(y+\beta) \), for integer values and halves.

Determine the value of \(\Gamma(n)\) using:

\( \Gamma(n) = (n-1)! \mbox{ or } \Gamma(n+1/2) = (n-1/2)\Gamma(n-1/2)\)

where \( \Gamma(1/2) = \sqrt(\pi)\)

Attempts simplification in cases like \( \Gamma(x+1)/\Gamma(x-1) = x*(x-1) \)

Definition at line 1371 of file Polylib.cpp.

1373 {
1374  double gamma = 1.0;
1375  double halfa = fabs(alpha - int(alpha));
1376  double halfb = fabs(beta - int(beta));
1377  if (halfa == 0.0 && halfb == 0.0) // integer value
1378  {
1379  int X = x + alpha;
1380  int Y = y + beta;
1381  if (X > Y)
1382  {
1383  for (int tmp = X - 1; tmp > Y - 1; tmp -= 1)
1384  gamma *= tmp;
1385  }
1386  else if (Y > X)
1387  {
1388  for (int tmp = Y - 1; tmp > X - 1; tmp -= 1)
1389  gamma *= tmp;
1390  gamma = 1. / gamma;
1391  }
1392  }
1393  else if (halfa == 0.5 && halfb == 0.5) // both are halves
1394  {
1395  double X = x + alpha;
1396  double Y = y + beta;
1397  if (X > Y)
1398  {
1399  for (int tmp = int(X); tmp > int(Y); tmp -= 1)
1400  gamma *= tmp - 0.5;
1401  }
1402  else if (Y > X)
1403  {
1404  for (int tmp = int(Y); tmp > int(X); tmp -= 1)
1405  gamma *= tmp - 0.5;
1406  gamma = 1. / gamma;
1407  }
1408  }
1409  else
1410  {
1411  double X = x + alpha;
1412  double Y = y + beta;
1413  while (X > 1 || Y > 1)
1414  {
1415  if (X > 1)
1416  {
1417  gamma *= X - 1.;
1418  X -= 1.;
1419  }
1420  if (Y > 1)
1421  {
1422  gamma /= Y - 1.;
1423  Y -= 1.;
1424  }
1425  }
1426  if (X == 0.5)
1427  {
1428  gamma *= sqrt(M_PI);
1429  }
1430  else if (Y == 0.5)
1431  {
1432  gamma /= sqrt(M_PI);
1433  }
1434  else
1435  {
1436  fprintf(stderr, "%lf or %lf is not of integer or half order\n", X,
1437  Y);
1438  }
1439  }
1440 
1441  return gamma;
1442 }

References Nektar::LibUtilities::beta, and tinysimd::sqrt().

Referenced by Dglj(), Dgrjm(), Dgrjp(), test_gamma_fraction(), zwgj(), zwglj(), zwgrjm(), and zwgrjp().

◆ hgj()

double Polylib::hgj ( const int  i,
const double  z,
const double *  zgj,
const int  np,
const double  alpha,
const double  beta 
)

Compute the value of the i th Lagrangian interpolant through the np Gauss-Jacobi points zgj at the arbitrary location z.

  • \( -1 \leq z \leq 1 \)
  • Uses the defintion of the Lagrangian interpolant:
    % \( \begin{array}{rcl} h_j(z) = \left\{ \begin{array}{ll} \displaystyle \frac{P_{np}^{\alpha,\beta}(z)} {[P_{np}^{\alpha,\beta}(z_j)]^\prime (z-z_j)} & \mbox{if $z \ne z_j$}\\ & \\ 1 & \mbox{if $z=z_j$} \end{array} \right. \end{array} \)

Definition at line 858 of file Polylib.cpp.

860 {
861  boost::ignore_unused(alpha, beta);
862  double zi, dz;
863 
864  zi = *(zgj + i);
865  dz = z - zi;
866  if (fabs(dz) < EPS)
867  return 1.0;
868 
869  return laginterp(z, i, zgj, np);
870 }
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:47
Polylib::laginterp
double laginterp(double z, int j, const double *zj, int np)
Definition: Polylib.cpp:87

References Nektar::LibUtilities::beta, EPS, and laginterp().

Referenced by Nektar::LibUtilities::Basis::GenBasis(), and Imgj().

◆ hglj()

double Polylib::hglj ( const int  i,
const double  z,
const double *  zglj,
const int  np,
const double  alpha,
const double  beta 
)

Compute the value of the i th Lagrangian interpolant through the np Gauss-Lobatto-Jacobi points zgrj at the arbitrary location z.

  • \( -1 \leq z \leq 1 \)
  • Uses the defintion of the Lagrangian interpolant:
    % \( \begin{array}{rcl} h_j(z) = \left\{ \begin{array}{ll} \displaystyle \frac{(1-z^2) P_{np-2}^{\alpha+1,\beta+1}(z)} {((1-z^2_j) [P_{np-2}^{\alpha+1,\beta+1}(z_j)]^\prime - 2 z_j P_{np-2}^{\alpha+1,\beta+1}(z_j) ) (z-z_j)}&\mbox{if $z \ne z_j$}\\ & \\ 1 & \mbox{if $z=z_j$} \end{array} \right. \end{array} \)

Definition at line 964 of file Polylib.cpp.

966 {
967  boost::ignore_unused(alpha, beta);
968 
969  double zi, dz;
970 
971  zi = *(zglj + i);
972  dz = z - zi;
973  if (fabs(dz) < EPS)
974  return 1.0;
975 
976  return laginterp(z, i, zglj, np);
977 }

References Nektar::LibUtilities::beta, EPS, and laginterp().

Referenced by Nektar::LibUtilities::Basis::GenBasis(), and Imglj().

◆ hgrjm()

double Polylib::hgrjm ( const int  i,
const double  z,
const double *  zgrj,
const int  np,
const double  alpha,
const double  beta 
)

Compute the value of the i th Lagrangian interpolant through the np Gauss-Radau-Jacobi points zgrj at the arbitrary location z. This routine assumes zgrj includes the point -1.

  • \( -1 \leq z \leq 1 \)
  • Uses the defintion of the Lagrangian interpolant:
    % \( \begin{array}{rcl} h_j(z) = \left\{ \begin{array}{ll} \displaystyle \frac{(1+z) P_{np-1}^{\alpha,\beta+1}(z)} {((1+z_j) [P_{np-1}^{\alpha,\beta+1}(z_j)]^\prime + P_{np-1}^{\alpha,\beta+1}(z_j) ) (z-z_j)} & \mbox{if $z \ne z_j$}\\ & \\ 1 & \mbox{if $z=z_j$} \end{array} \right. \end{array} \)

Definition at line 893 of file Polylib.cpp.

895 {
896  boost::ignore_unused(alpha, beta);
897 
898  double zi, dz;
899 
900  zi = *(zgrj + i);
901  dz = z - zi;
902  if (fabs(dz) < EPS)
903  return 1.0;
904 
905  return laginterp(z, i, zgrj, np);
906 }

References Nektar::LibUtilities::beta, EPS, and laginterp().

Referenced by Imgrjm().

◆ hgrjp()

double Polylib::hgrjp ( const int  i,
const double  z,
const double *  zgrj,
const int  np,
const double  alpha,
const double  beta 
)

Compute the value of the i th Lagrangian interpolant through the np Gauss-Radau-Jacobi points zgrj at the arbitrary location z. This routine assumes zgrj includes the point +1.

  • \( -1 \leq z \leq 1 \)
  • Uses the defintion of the Lagrangian interpolant:
    % \( \begin{array}{rcl} h_j(z) = \left\{ \begin{array}{ll} \displaystyle \frac{(1-z) P_{np-1}^{\alpha+1,\beta}(z)} {((1-z_j) [P_{np-1}^{\alpha+1,\beta}(z_j)]^\prime - P_{np-1}^{\alpha+1,\beta}(z_j) ) (z-z_j)} & \mbox{if $z \ne z_j$}\\ & \\ 1 & \mbox{if $z=z_j$} \end{array} \right. \end{array} \)

Definition at line 929 of file Polylib.cpp.

931 {
932  boost::ignore_unused(alpha, beta);
933  double zi, dz;
934 
935  zi = *(zgrj + i);
936  dz = z - zi;
937  if (fabs(dz) < EPS)
938  return 1.0;
939 
940  return laginterp(z, i, zgrj, np);
941 }

References Nektar::LibUtilities::beta, EPS, and laginterp().

Referenced by Imgrjp().

◆ ImagBesselComp()

std::complex< Nektar::NekDouble > Polylib::ImagBesselComp ( int  n,
std::complex< Nektar::NekDouble >  y 
)

Calcualte the bessel function of the first kind with complex double input y. Taken from Numerical Recipies in C.

Returns a complex double

Definition at line 1837 of file Polylib.cpp.

1839 {
1840  std::complex<Nektar::NekDouble> z(1.0, 0.0);
1841  std::complex<Nektar::NekDouble> zbes(1.0, 0.0);
1842  std::complex<Nektar::NekDouble> zarg;
1843  Nektar::NekDouble tol = 1e-15;
1844  int maxit = 10000;
1845  int i = 1;
1846 
1847  zarg = -0.25 * y * y;
1848 
1849  while (abs(z) > tol && i <= maxit)
1850  {
1851  z = z * (1.0 / i / (i + n) * zarg);
1852  if (abs(z) <= tol)
1853  break;
1854  zbes = zbes + z;
1855  i++;
1856  }
1857  zarg = 0.5 * y;
1858  for (i = 1; i <= n; i++)
1859  {
1860  zbes = zbes * zarg;
1861  }
1862  return zbes;
1863 }
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
tinysimd::abs
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:298

References tinysimd::abs().

Referenced by Nektar::IncNavierStokes::SetUpWomersley().

◆ Imgj()

void Polylib::Imgj ( double *  im,
const double *  zgj,
const double *  zm,
const int  nz,
const int  mz,
const double  alpha,
const double  beta 
)

Interpolation Operator from Gauss-Jacobi points to an arbitrary distribution at points zm.

  • Computes the one-dimensional interpolation matrix, im, to interpolate a function from at Gauss-Jacobi distribution of nz zeros zgrj to an arbitrary distribution of mz points zm, i.e.
    \( u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \)

Definition at line 992 of file Polylib.cpp.

994 {
995  double zp;
996  int i, j;
997 
998  for (i = 0; i < nz; ++i)
999  {
1000  for (j = 0; j < mz; ++j)
1001  {
1002  zp = zm[j];
1003  im[i * mz + j] = hgj(i, zp, zgj, nz, alpha, beta);
1004  }
1005  }
1006 
1007  return;
1008 }
Polylib::hgj
double hgj(const int i, const double z, const double *zgj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the np Gauss-Jacobi points zgj at the ar...
Definition: Polylib.cpp:858

References Nektar::LibUtilities::beta, and hgj().

Referenced by Nektar::LibUtilities::GaussPoints::CalculateInterpMatrix(), and main().

◆ Imglj()

void Polylib::Imglj ( double *  im,
const double *  zglj,
const double *  zm,
const int  nz,
const int  mz,
const double  alpha,
const double  beta 
)

Interpolation Operator from Gauss-Lobatto-Jacobi points to an arbitrary distrubtion at points zm.

  • Computes the one-dimensional interpolation matrix, im, to interpolate a function from at Gauss-Lobatto-Jacobi distribution of nz zeros zgrj (where zgrj[0]=-1) to an arbitrary distribution of mz points zm, i.e.
    \( u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \)

Definition at line 1085 of file Polylib.cpp.

1087 {
1088  double zp;
1089  int i, j;
1090 
1091  for (i = 0; i < nz; i++)
1092  {
1093  for (j = 0; j < mz; j++)
1094  {
1095  zp = zm[j];
1096  im[i * mz + j] = hglj(i, zp, zglj, nz, alpha, beta);
1097  }
1098  }
1099 
1100  return;
1101 }
Polylib::hglj
double hglj(const int i, const double z, const double *zglj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the np Gauss-Lobatto-Jacobi points zgrj ...
Definition: Polylib.cpp:964

References Nektar::LibUtilities::beta, and hglj().

Referenced by Nektar::LibUtilities::GaussPoints::CalculateInterpMatrix(), and main().

◆ Imgrjm()

void Polylib::Imgrjm ( double *  im,
const double *  zgrj,
const double *  zm,
const int  nz,
const int  mz,
const double  alpha,
const double  beta 
)

Interpolation Operator from Gauss-Radau-Jacobi points (including z=-1) to an arbitrary distrubtion at points zm.

  • Computes the one-dimensional interpolation matrix, im, to interpolate a function from at Gauss-Radau-Jacobi distribution of nz zeros zgrj (where zgrj[0]=-1) to an arbitrary distribution of mz points zm, i.e.
    \( u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \)

Definition at line 1023 of file Polylib.cpp.

1025 {
1026  double zp;
1027  int i, j;
1028 
1029  for (i = 0; i < nz; i++)
1030  {
1031  for (j = 0; j < mz; j++)
1032  {
1033  zp = zm[j];
1034  im[i * mz + j] = hgrjm(i, zp, zgrj, nz, alpha, beta);
1035  }
1036  }
1037 
1038  return;
1039 }
Polylib::hgrjm
double hgrjm(const int i, const double z, const double *zgrj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the np Gauss-Radau-Jacobi points zgrj at...
Definition: Polylib.cpp:893

References Nektar::LibUtilities::beta, and hgrjm().

Referenced by Nektar::LibUtilities::GaussPoints::CalculateInterpMatrix(), and main().

◆ Imgrjp()

void Polylib::Imgrjp ( double *  im,
const double *  zgrj,
const double *  zm,
const int  nz,
const int  mz,
const double  alpha,
const double  beta 
)

Interpolation Operator from Gauss-Radau-Jacobi points (including z=1) to an arbitrary distrubtion at points zm.

  • Computes the one-dimensional interpolation matrix, im, to interpolate a function from at Gauss-Radau-Jacobi distribution of nz zeros zgrj (where zgrj[nz-1]=1) to an arbitrary distribution of mz points zm, i.e.
    \( u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \)

Definition at line 1054 of file Polylib.cpp.

1056 {
1057  double zp;
1058  int i, j;
1059 
1060  for (i = 0; i < nz; i++)
1061  {
1062  for (j = 0; j < mz; j++)
1063  {
1064  zp = zm[j];
1065  im[i * mz + j] = hgrjp(i, zp, zgrj, nz, alpha, beta);
1066  }
1067  }
1068 
1069  return;
1070 }
Polylib::hgrjp
double hgrjp(const int i, const double z, const double *zgrj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the np Gauss-Radau-Jacobi points zgrj at...
Definition: Polylib.cpp:929

References Nektar::LibUtilities::beta, and hgrjp().

Referenced by Nektar::LibUtilities::GaussPoints::CalculateInterpMatrix(), and main().

◆ jacobd()

void Polylib::jacobd ( const int  np,
const double *  z,
double *  polyd,
const int  n,
const double  alpha,
const double  beta 
)

Calculate the derivative of Jacobi polynomials.

  • Generates a vector poly of values of the derivative of the n th order Jacobi polynomial \( P^(\alpha,\beta)_n(z)\) at the np points z.
  • To do this we have used the relation
    \( \frac{d}{dz} P^{\alpha,\beta}_n(z) = \frac{1}{2} (\alpha + \beta + n + 1) P^{\alpha,\beta}_n(z) \)
  • This formulation is valid for \( -1 \leq z \leq 1 \)

Definition at line 1293 of file Polylib.cpp.

1295 {
1296  int i;
1297  double one = 1.0;
1298  if (n == 0)
1299  for (i = 0; i < np; ++i)
1300  polyd[i] = 0.0;
1301  else
1302  {
1303  // jacobf(np,z,polyd,n-1,alpha+one,beta+one);
1304  jacobfd(np, z, polyd, NULL, n - 1, alpha + one, beta + one);
1305  for (i = 0; i < np; ++i)
1306  polyd[i] *= 0.5 * (alpha + beta + (double)n + one);
1307  }
1308  return;
1309 }
Polylib::jacobfd
void jacobfd(const int np, const double *z, double *poly_in, double *polyd, const int n, const double alpha, const double beta)
Routine to calculate Jacobi polynomials, , and their first derivative, .
Definition: Polylib.cpp:1181

References Nektar::LibUtilities::beta, and jacobfd().

Referenced by Dgj(), Dglj(), Dgrjm(), Dgrjp(), Nektar::SolverUtils::AdvectionFR::SetupCFunctions(), Nektar::SolverUtils::DiffusionLFR::SetupCFunctions(), Nektar::SolverUtils::DiffusionLFRNS::SetupCFunctions(), Nektar::LibUtilities::NodalUtilTriangle::v_OrthoBasisDeriv(), Nektar::LibUtilities::NodalUtilTetrahedron::v_OrthoBasisDeriv(), Nektar::LibUtilities::NodalUtilPrism::v_OrthoBasisDeriv(), Nektar::LibUtilities::NodalUtilQuad::v_OrthoBasisDeriv(), Nektar::LibUtilities::NodalUtilHex::v_OrthoBasisDeriv(), and zwgj().

◆ jacobfd()

void Polylib::jacobfd ( const int  np,
const double *  z,
double *  poly_in,
double *  polyd,
const int  n,
const double  alpha,
const double  beta 
)

Routine to calculate Jacobi polynomials, \( P^{\alpha,\beta}_n(z) \), and their first derivative, \( \frac{d}{dz} P^{\alpha,\beta}_n(z) \).

  • This function returns the vectors poly_in and poly_d containing the value of the \( n^th \) order Jacobi polynomial \( P^{\alpha,\beta}_n(z) \alpha > -1, \beta > -1 \) and its derivative at the np points in z[i]
  • If poly_in = NULL then only calculate derivatice
  • If polyd = NULL then only calculate polynomial
  • To calculate the polynomial this routine uses the recursion relationship (see appendix A ref [4]) : \( \begin{array}{rcl} P^{\alpha,\beta}_0(z) &=& 1 \\ P^{\alpha,\beta}_1(z) &=& \frac{1}{2} [ \alpha-\beta+(\alpha+\beta+2)z] \\ a^1_n P^{\alpha,\beta}_{n+1}(z) &=& (a^2_n + a^3_n z) P^{\alpha,\beta}_n(z) - a^4_n P^{\alpha,\beta}_{n-1}(z) \\ a^1_n &=& 2(n+1)(n+\alpha + \beta + 1)(2n + \alpha + \beta) \\ a^2_n &=& (2n + \alpha + \beta + 1)(\alpha^2 - \beta^2) \\ a^3_n &=& (2n + \alpha + \beta)(2n + \alpha + \beta + 1) (2n + \alpha + \beta + 2) \\ a^4_n &=& 2(n+\alpha)(n+\beta)(2n + \alpha + \beta + 2) \end{array} \)
  • To calculate the derivative of the polynomial this routine uses the relationship (see appendix A ref [4]) : \( \begin{array}{rcl} b^1_n(z)\frac{d}{dz} P^{\alpha,\beta}_n(z)&=&b^2_n(z)P^{\alpha,\beta}_n(z) + b^3_n(z) P^{\alpha,\beta}_{n-1}(z) \hspace{2.2cm} \\ b^1_n(z) &=& (2n+\alpha + \beta)(1-z^2) \\ b^2_n(z) &=& n[\alpha - \beta - (2n+\alpha + \beta)z]\\ b^3_n(z) &=& 2(n+\alpha)(n+\beta) \end{array} \)
  • Note the derivative from this routine is only valid for -1 < z < 1.

Definition at line 1181 of file Polylib.cpp.

1183 {
1184  int i;
1185  double zero = 0.0, one = 1.0, two = 2.0;
1186 
1187  if (!np)
1188  return;
1189 
1190  if (n == 0)
1191  {
1192  if (poly_in)
1193  for (i = 0; i < np; ++i)
1194  poly_in[i] = one;
1195  if (polyd)
1196  for (i = 0; i < np; ++i)
1197  polyd[i] = zero;
1198  }
1199  else if (n == 1)
1200  {
1201  if (poly_in)
1202  for (i = 0; i < np; ++i)
1203  poly_in[i] = 0.5 * (alpha - beta + (alpha + beta + two) * z[i]);
1204  if (polyd)
1205  for (i = 0; i < np; ++i)
1206  polyd[i] = 0.5 * (alpha + beta + two);
1207  }
1208  else
1209  {
1210  int k;
1211  double a1, a2, a3, a4;
1212  double two = 2.0, apb = alpha + beta;
1213  double *poly, *polyn1, *polyn2;
1214 
1215  if (poly_in)
1216  { // switch for case of no poynomial function return
1217  polyn1 = (double *)malloc(2 * np * sizeof(double));
1218  polyn2 = polyn1 + np;
1219  poly = poly_in;
1220  }
1221  else
1222  {
1223  polyn1 = (double *)malloc(3 * np * sizeof(double));
1224  polyn2 = polyn1 + np;
1225  poly = polyn2 + np;
1226  }
1227 
1228  for (i = 0; i < np; ++i)
1229  {
1230  polyn2[i] = one;
1231  polyn1[i] = 0.5 * (alpha - beta + (alpha + beta + two) * z[i]);
1232  }
1233 
1234  for (k = 2; k <= n; ++k)
1235  {
1236  a1 = two * k * (k + apb) * (two * k + apb - two);
1237  a2 = (two * k + apb - one) * (alpha * alpha - beta * beta);
1238  a3 =
1239  (two * k + apb - two) * (two * k + apb - one) * (two * k + apb);
1240  a4 = two * (k + alpha - one) * (k + beta - one) * (two * k + apb);
1241 
1242  a2 /= a1;
1243  a3 /= a1;
1244  a4 /= a1;
1245 
1246  for (i = 0; i < np; ++i)
1247  {
1248  poly[i] = (a2 + a3 * z[i]) * polyn1[i] - a4 * polyn2[i];
1249  polyn2[i] = polyn1[i];
1250  polyn1[i] = poly[i];
1251  }
1252  }
1253 
1254  if (polyd)
1255  {
1256  a1 = n * (alpha - beta);
1257  a2 = n * (two * n + alpha + beta);
1258  a3 = two * (n + alpha) * (n + beta);
1259  a4 = (two * n + alpha + beta);
1260  a1 /= a4;
1261  a2 /= a4;
1262  a3 /= a4;
1263 
1264  // note polyn2 points to polyn1 at end of poly iterations
1265  for (i = 0; i < np; ++i)
1266  {
1267  polyd[i] = (a1 - a2 * z[i]) * poly[i] + a3 * polyn2[i];
1268  polyd[i] /= (one - z[i] * z[i]);
1269  }
1270  }
1271 
1272  free(polyn1);
1273  }
1274 
1275  return;
1276 }

References Nektar::LibUtilities::beta.

Referenced by Nektar::LibUtilities::Basis::GenBasis(), jacobd(), Jacobz(), main(), Nektar::LibUtilities::NodalUtilTriangle::v_OrthoBasis(), Nektar::LibUtilities::NodalUtilTetrahedron::v_OrthoBasis(), Nektar::LibUtilities::NodalUtilPrism::v_OrthoBasis(), Nektar::LibUtilities::NodalUtilQuad::v_OrthoBasis(), Nektar::LibUtilities::NodalUtilHex::v_OrthoBasis(), Nektar::LibUtilities::NodalUtilTriangle::v_OrthoBasisDeriv(), Nektar::LibUtilities::NodalUtilTetrahedron::v_OrthoBasisDeriv(), Nektar::LibUtilities::NodalUtilPrism::v_OrthoBasisDeriv(), Nektar::LibUtilities::NodalUtilQuad::v_OrthoBasisDeriv(), Nektar::LibUtilities::NodalUtilHex::v_OrthoBasisDeriv(), zwglj(), zwgrjm(), and zwgrjp().

◆ Jacobz()

static void Polylib::Jacobz ( const int  n,
double *  z,
const double  alpha,
const double  beta 
)
static

Calculate the n zeros, z, of the Jacobi polynomial, i.e. \( P_n^{\alpha,\beta}(z) = 0 \).

This routine is only value for \(( \alpha > -1, \beta > -1)\) and uses polynomial deflation in a Newton iteration

Definition at line 1452 of file Polylib.cpp.

1454 {
1455  int i, j, k;
1456  double dth = M_PI / (2.0 * (double)n);
1457  double poly, pder, rlast = 0.0;
1458  double sum, delr, r;
1459  double one = 1.0, two = 2.0;
1460 
1461  if (!n)
1462  return;
1463 
1464  for (k = 0; k < n; ++k)
1465  {
1466  r = -cos((two * (double)k + one) * dth);
1467  if (k)
1468  r = 0.5 * (r + rlast);
1469 
1470  for (j = 1; j < STOP; ++j)
1471  {
1472  jacobfd(1, &r, &poly, &pder, n, alpha, beta);
1473 
1474  for (i = 0, sum = 0.0; i < k; ++i)
1475  sum += one / (r - z[i]);
1476 
1477  delr = -poly / (pder - sum * poly);
1478  r += delr;
1479  if (fabs(delr) < EPS)
1480  break;
1481  }
1482  z[k] = r;
1483  rlast = r;
1484  }
1485  return;
1486 }
STOP
#define STOP
Maximum number of iterations in polynomial defalation routine Jacobz.
Definition: Polylib.cpp:45

References Nektar::LibUtilities::beta, EPS, jacobfd(), and STOP.

◆ JacZeros()

void Polylib::JacZeros ( const int  n,
double *  a,
double *  b,
const double  alpha,
const double  beta 
)

Zero and Weight determination through the eigenvalues and eigenvectors of a tridiagonal matrix from the three term recurrence relationship.

Set up a symmetric tridiagonal matrix

\( \left [ \begin{array}{ccccc} a[0] & b[0] & & & \\ b[0] & a[1] & b[1] & & \\ 0 & \ddots & \ddots & \ddots & \\ & & \ddots & \ddots & b[n-2] \\ & & & b[n-2] & a[n-1] \end{array} \right ] \)

Where the coefficients a[n], b[n] come from the recurrence relation

\( b_j p_j(z) = (z - a_j ) p_{j-1}(z) - b_{j-1} p_{j-2}(z) \)

where \( j=n+1\) and \(p_j(z)\) are the Jacobi (normalized) orthogonal polynomials \( \alpha,\beta > -1\)( integer values and halves). Since the polynomials are orthonormalized, the tridiagonal matrix is guaranteed to be symmetric. The eigenvalues of this matrix are the zeros of the Jacobi polynomial.

Definition at line 1512 of file Polylib.cpp.

1514 {
1515 
1516  int i, j;
1517  RecCoeff(n, a, b, alpha, beta);
1518 
1519  double **z = new double *[n];
1520  for (i = 0; i < n; i++)
1521  {
1522  z[i] = new double[n];
1523  for (j = 0; j < n; j++)
1524  {
1525  z[i][j] = 0.0;
1526  }
1527  }
1528  for (i = 0; i < n; i++)
1529  {
1530  z[i][i] = 1.0;
1531  }
1532 
1533  // find eigenvalues and eigenvectors
1534  TriQL(n, a, b, z);
1535 
1536  delete[] z;
1537  return;
1538 }
Polylib::RecCoeff
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and b of the orthogonal polynomials.
Definition: Polylib.cpp:1544
Polylib::TriQL
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
Definition: Polylib.cpp:1603

References Nektar::LibUtilities::beta, RecCoeff(), and TriQL().

◆ JKMatrix()

void Polylib::JKMatrix ( int  n,
double *  a,
double *  b 
)

Calcualtes the Jacobi-kronrod matrix by determining the a and coefficients.

The first 3n+1 coefficients are already known

For more information refer to: "Dirk P. Laurie, Calcualtion of Gauss-Kronrod quadrature rules"

Definition at line 1716 of file Polylib.cpp.

1717 {
1718  int i, j, k, m;
1719  // Working storage
1720  int size = (int)floor(n / 2.0) + 2;
1721  double *s = new double[size];
1722  double *t = new double[size];
1723 
1724  // Initialize s and t to zero
1725  for (i = 0; i < size; i++)
1726  {
1727  s[i] = 0.0;
1728  t[i] = 0.0;
1729  }
1730 
1731  t[1] = b[n + 1];
1732  for (m = 0; m <= n - 2; m++)
1733  {
1734  double u = 0.0;
1735  for (k = (int)floor((m + 1) / 2.0); k >= 0; k--)
1736  {
1737  int l = m - k;
1738  u = u + (a[k + n + 1] - a[l]) * t[k + 1] + b[k + n + 1] * s[k] -
1739  b[l] * s[k + 1];
1740  s[k + 1] = u;
1741  }
1742 
1743  // Swap the contents of s and t
1744  double *hold = s;
1745  s = t;
1746  t = hold;
1747  }
1748 
1749  for (j = (int)floor(n / 2.0); j >= 0; j--)
1750  {
1751  s[j + 1] = s[j];
1752  }
1753 
1754  for (m = n - 1; m <= 2 * n - 3; m++)
1755  {
1756  double u = 0;
1757  for (k = m + 1 - n; k <= floor((m - 1) / 2.0); k++)
1758  {
1759  int l = m - k;
1760  j = n - 1 - l;
1761  u = u - (a[k + n + 1] - a[l]) * t[j + 1] - b[k + n + 1] * s[j + 1] +
1762  b[l] * s[j + 2];
1763  s[j + 1] = u;
1764  }
1765 
1766  if (m % 2 == 0)
1767  {
1768  k = m / 2;
1769  a[k + n + 1] =
1770  a[k] + (s[j + 1] - b[k + n + 1] * s[j + 2]) / t[j + 2];
1771  }
1772  else
1773  {
1774  k = (m + 1) / 2;
1775  b[k + n + 1] = s[j + 1] / s[j + 2];
1776  }
1777 
1778  // Swap the contents of s and t
1779  double *hold = s;
1780  s = t;
1781  t = hold;
1782  }
1783 
1784  a[2 * n] = a[n - 1] - b[2 * n] * s[1] / t[1];
1785 }

Referenced by zwgk(), zwlk(), and zwrk().

◆ laginterp()

double Polylib::laginterp ( double  z,
int  j,
const double *  zj,
int  np 
)

Definition at line 87 of file Polylib.cpp.

88 {
89  double temp = 1.0;
90  for (int i = 0; i < np; i++)
91  {
92  if (j != i)
93  {
94  temp *= optdiff(z, zj[i]) / (zj[j] - zj[i]);
95  }
96  }
97  return temp;
98 }
Polylib::optdiff
double optdiff(double xl, double xr)
The following function is used to circumvent/reduce "Subtractive Cancellation" The expression 1/dz is...
Definition: Polylib.cpp:57

References optdiff().

Referenced by hgj(), hglj(), hgrjm(), and hgrjp().

◆ optdiff()

double Polylib::optdiff ( double  xl,
double  xr 
)

The following function is used to circumvent/reduce "Subtractive Cancellation" The expression 1/dz is replaced by optinvsub(.,.) Added on 26 April 2017.

Definition at line 57 of file Polylib.cpp.

58 {
59  double m_xln, m_xrn;
60  int m_expn;
61  int m_digits = static_cast<int>(fabs(floor(log10(DBL_EPSILON))) - 1);
62 
63  if (fabs(xl - xr) < 1.e-4)
64  {
65 
66  m_expn = static_cast<int>(floor(log10(fabs(xl - xr))));
67  m_xln = xl * powl(10.0L, -m_expn) -
68  floor(xl * powl(10.0L,
69  -m_expn)); // substract the digits overlap part
70  m_xrn =
71  xr * powl(10.0L, -m_expn) -
72  floor(xl *
73  powl(10.0L,
74  -m_expn)); // substract the common digits overlap part
75  m_xln =
76  round(m_xln * powl(10.0L, m_digits + m_expn)); // git rid of rubbish
77  m_xrn = round(m_xrn * powl(10.0L, m_digits + m_expn));
78 
79  return powl(10.0L, -m_digits) * (m_xln - m_xrn);
80  }
81  else
82  {
83  return (xl - xr);
84  }
85 }

Referenced by laginterp().

◆ polycoeffs()

void Polylib::polycoeffs ( const int  i,
const double *  z,
double *  c,
const int  np 
)

Definition at line 1103 of file Polylib.cpp.

1104 {
1105  int j, k, m;
1106 
1107  // Compute denominator
1108  double d = 1.0;
1109  for (j = 0; j < np; j++)
1110  {
1111  if (i != j)
1112  {
1113  d *= z[i] - z[j];
1114  }
1115  c[j] = 0.0;
1116  }
1117 
1118  // Compute coefficient
1119  c[0] = 1.0;
1120  m = 0;
1121  for (j = 0; j < np; j++)
1122  {
1123  if (i != j)
1124  {
1125  m += 1;
1126  c[m] = c[m - 1];
1127  for (k = m - 1; k > 0; k--)
1128  {
1129  c[k] *= -1.0 * z[j];
1130  c[k] += c[k - 1];
1131  }
1132  c[0] *= -1.0 * z[j];
1133  }
1134  }
1135  for (j = 0; j < np; j++)
1136  {
1137  c[j] /= d;
1138  }
1139 }

Referenced by Qg().

◆ Qg()

void Polylib::Qg ( double *  Q,
const double *  z,
const int  np,
const int  offset 
)

Compute the Integration Matrix.

Definition at line 584 of file Polylib.cpp.

585 {
586 
587  if (np <= 0)
588  {
589  Q[0] = 0.0;
590  }
591  else
592  {
593  int i, j, k;
594  double *pd;
595 
596  pd = (double *)malloc(np * sizeof(double));
597 
598  for (i = 0; i < np; i++)
599  {
600  polycoeffs(i, z, pd, np);
601  for (j = 0; j < np; j++)
602  {
603  Q[j * (np + offset) + i + offset] = 0.0;
604  for (k = 0; k < np; k++)
605  {
606  Q[j * (np + offset) + i + offset] +=
607  pd[k] / (k + 1) * pow(z[j], k + 1);
608  }
609  }
610  }
611  free(pd);
612  }
613  return;
614 }
Polylib::polycoeffs
void polycoeffs(const int i, const double *z, double *c, const int np)
Definition: Polylib.cpp:1103

References polycoeffs().

Referenced by main(), and Nektar::LibUtilities::TimeIntegrationSchemeSDC::v_InitializeScheme().

◆ RecCoeff()

static void Polylib::RecCoeff ( const int  n,
double *  a,
double *  b,
const double  alpha,
const double  beta 
)
static

The routine finds the recurrence coefficients a and b of the orthogonal polynomials.

Definition at line 1544 of file Polylib.cpp.

1546 {
1547 
1548  int i;
1549  double apb, apbi, a2b2;
1550 
1551  if (!n)
1552  return;
1553 
1554  // generate normalised terms
1555  apb = alpha + beta;
1556  apbi = 2.0 + apb;
1557 
1558  b[0] = pow(2.0, apb + 1.0) * gammaF(alpha + 1.0) * gammaF(beta + 1.0) /
1559  gammaF(apbi); // MuZero
1560  a[0] = (beta - alpha) / apbi;
1561  b[1] = (4.0 * (1.0 + alpha) * (1.0 + beta) / ((apbi + 1.0) * apbi * apbi));
1562 
1563  a2b2 = beta * beta - alpha * alpha;
1564 
1565  for (i = 1; i < n - 1; i++)
1566  {
1567  apbi = 2.0 * (i + 1) + apb;
1568  a[i] = a2b2 / ((apbi - 2.0) * apbi);
1569  b[i + 1] = (4.0 * (i + 1) * (i + 1 + alpha) * (i + 1 + beta) *
1570  (i + 1 + apb) / ((apbi * apbi - 1) * apbi * apbi));
1571  }
1572 
1573  apbi = 2.0 * n + apb;
1574  a[n - 1] = a2b2 / ((apbi - 2.0) * apbi);
1575 }

References Nektar::LibUtilities::beta, and gammaF().

Referenced by JacZeros(), zwgk(), zwlk(), and zwrk().

◆ TriQL()

static void Polylib::TriQL ( const int  n,
double *  d,
double *  e,
double **  z 
)
static

QL algorithm for symmetric tridiagonal matrix.

This subroutine is a translation of an algol procedure, num. math. 12, 377-383(1968) by martin and wilkinson, as modified in num. math. 15, 450(1970) by dubrulle. Handbook for auto. comp., vol.ii-linear algebra, 241-248(1971). This is a modified version from numerical recipes.

This subroutine finds the eigenvalues and first components of the eigenvectors of a symmetric tridiagonal matrix by the implicit QL method.

on input:

  • n is the order of the matrix;
  • d contains the diagonal elements of the input matrix;
  • e contains the subdiagonal elements of the input matrix in its first n-2 positions.
    • z is the n by n identity matrix

on output:

  • d contains the eigenvalues in ascending order.
  • e contains the weight values - modifications of the first component of normalised eigenvectors

Definition at line 1603 of file Polylib.cpp.

1604 {
1605  int m, l, iter, i, k;
1606  double s, r, p, g, f, dd, c, b;
1607 
1608  double MuZero = e[0];
1609 
1610  // Renumber the elements of e
1611  for (i = 0; i < n - 1; i++)
1612  {
1613  e[i] = sqrt(e[i + 1]);
1614  }
1615  e[n - 1] = 0.0;
1616 
1617  for (l = 0; l < n; l++)
1618  {
1619  iter = 0;
1620  do
1621  {
1622  for (m = l; m < n - 1; m++)
1623  {
1624  dd = fabs(d[m]) + fabs(d[m + 1]);
1625  if (fabs(e[m]) + dd == dd)
1626  break;
1627  }
1628  if (m != l)
1629  {
1630  if (iter++ == STOP)
1631  {
1632  fprintf(stderr, "triQL: Too many iterations in TQLI");
1633  exit(1);
1634  }
1635  g = (d[l + 1] - d[l]) / (2.0 * e[l]);
1636  r = sqrt((g * g) + 1.0);
1637  g = d[m] - d[l] + e[l] / (g + sign(r, g));
1638  s = c = 1.0;
1639  p = 0.0;
1640  for (i = m - 1; i >= l; i--)
1641  {
1642  f = s * e[i];
1643  b = c * e[i];
1644  if (fabs(f) >= fabs(g))
1645  {
1646  c = g / f;
1647  r = sqrt((c * c) + 1.0);
1648  e[i + 1] = f * r;
1649  c *= (s = 1.0 / r);
1650  }
1651  else
1652  {
1653  s = f / g;
1654  r = sqrt((s * s) + 1.0);
1655  e[i + 1] = g * r;
1656  s *= (c = 1.0 / r);
1657  }
1658  g = d[i + 1] - p;
1659  r = (d[i] - g) * s + 2.0 * c * b;
1660  p = s * r;
1661  d[i + 1] = g + p;
1662  g = c * r - b;
1663 
1664  // Calculate the eigenvectors
1665  for (k = 0; k < n; k++)
1666  {
1667  f = z[k][i + 1];
1668  z[k][i + 1] = s * z[k][i] + c * f;
1669  z[k][i] = c * z[k][i] - s * f;
1670  }
1671  }
1672  d[l] = d[l] - p;
1673  e[l] = g;
1674  e[m] = 0.0;
1675  }
1676  } while (m != l);
1677  }
1678 
1679  // order eigenvalues
1680  // Since we only need the first component of the eigenvectors
1681  // to calcualte the weight, we only swap the first components
1682  for (i = 0; i < n - 1; ++i)
1683  {
1684  k = i;
1685  p = d[i];
1686  for (l = i + 1; l < n; ++l)
1687  if (d[l] < p)
1688  {
1689  k = l;
1690  p = d[l];
1691  }
1692  d[k] = d[i];
1693  d[i] = p;
1694 
1695  double temp = z[0][k];
1696  z[0][k] = z[0][i];
1697  z[0][i] = temp;
1698  }
1699 
1700  // Calculate the weights
1701  for (i = 0; i < n; i++)
1702  {
1703  e[i] = MuZero * z[0][i] * z[0][i];
1704  }
1705 }
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:49
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350

References CellMLToNektar.cellml_metadata::p, sign, tinysimd::sqrt(), and STOP.

Referenced by JacZeros(), zwgk(), zwlk(), and zwrk().

◆ zwgj()

void Polylib::zwgj ( double *  z,
double *  w,
const int  np,
const double  alpha,
const double  beta 
)

Gauss-Jacobi zeros and weights.

  • Generate np Gauss Jacobi zeros, z, and weights,w, associated with the Jacobi polynomial \( P^{\alpha,\beta}_{np}(z) \),
  • Exact for polynomials of order 2np-1 or less

Definition at line 133 of file Polylib.cpp.

135 {
136  int i;
137  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
138 
139  jacobz(np, z, alpha, beta);
140  jacobd(np, z, w, np, alpha, beta);
141 
142  fac = pow(two, apb + one) * gammaFracGammaF(np + 1, alpha, np + 1, 0.0) *
143  gammaFracGammaF(np + 1, beta, np + 1, apb);
144 
145  for (i = 0; i < np; ++i)
146  w[i] = fac / (w[i] * w[i] * (one - z[i] * z[i]));
147 
148  return;
149 }
jacobz
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:104

References Nektar::LibUtilities::beta, gammaFracGammaF(), jacobd(), and jacobz.

Referenced by main(), and Nektar::LibUtilities::GaussPoints::v_CalculatePoints().

◆ zwgk()

void Polylib::zwgk ( double *  z,
double *  w,
const int  npt,
const double  alpha,
const double  beta 
)

Gauss-Kronrod-Jacobi zeros and weights.

  • Generate npt=2*np+1 Gauss-Kronrod Jacobi zeros, z, and weights,w, associated with the Jacobi polynomial \( P^{\alpha,\beta}_{np}(z) \),
  • Exact for polynomials of order 3np+1 or less

Definition at line 290 of file Polylib.cpp.

292 {
293 
294  int np = (npt - 1) / 2;
295 
296  int i, j;
297 
298  // number of kronrod points associated with the np gauss rule
299  int kpoints = 2 * np + 1;
300 
301  // Define the number of required recurrence coefficents
302  int ncoeffs = (int)floor(3.0 * (np + 1) / 2);
303 
304  // Define arrays for the recurrence coefficients
305  // We will use these arrays for the Kronrod results too, hence the
306  // reason for the size of the arrays
307  double *a = new double[kpoints];
308  double *b = new double[kpoints];
309 
310  // Initialize a and b to zero
311  for (i = 0; i < kpoints; i++)
312  {
313  a[i] = 0.0;
314  b[i] = 0.0;
315  }
316 
317  // Call the routine to calculate the recurrence coefficients
318  RecCoeff(ncoeffs, a, b, alpha, beta);
319 
320  // Call the routine to caluclate the jacobi-Kronrod matrix
321  JKMatrix(np, a, b);
322 
323  // Set up the identity matrix
324  double **zmatrix = new double *[kpoints];
325  for (i = 0; i < kpoints; i++)
326  {
327  zmatrix[i] = new double[kpoints];
328  for (j = 0; j < kpoints; j++)
329  {
330  zmatrix[i][j] = 0.0;
331  }
332  }
333  for (i = 0; i < kpoints; i++)
334  {
335  zmatrix[i][i] = 1.0;
336  }
337 
338  // Calculte the points and weights
339  TriQL(kpoints, a, b, zmatrix);
340 
341  for (i = 0; i < kpoints; i++)
342  {
343  z[i] = a[i];
344  w[i] = b[i];
345  }
346  delete[] a;
347  delete[] b;
348  for (i = 0; i < kpoints; i++)
349  {
350  delete[] zmatrix[i];
351  }
352  delete[] zmatrix;
353 }
Polylib::JKMatrix
void JKMatrix(int, double *, double *)
Calcualtes the Jacobi-kronrod matrix by determining the a and coefficients.
Definition: Polylib.cpp:1716

References Nektar::LibUtilities::beta, JKMatrix(), RecCoeff(), and TriQL().

Referenced by Nektar::LibUtilities::GaussPoints::v_CalculatePoints().

◆ zwglj()

void Polylib::zwglj ( double *  z,
double *  w,
const int  np,
const double  alpha,
const double  beta 
)

Gauss-Lobatto-Jacobi zeros and weights with end point at z=-1,1.

  • Generate np Gauss-Lobatto-Jacobi points, z, and weights, w, associated with polynomial \( (1-z)(1+z) P^{\alpha+1,\beta+1}_{np-2}(z) \)
  • Exact for polynomials of order 2np-3 or less

Definition at line 241 of file Polylib.cpp.

243 {
244 
245  if (np == 1)
246  {
247  z[0] = 0.0;
248  w[0] = 2.0;
249  }
250  else if (np == 2)
251  {
252  z[0] = -1.0;
253  z[1] = 1.0;
254 
255  w[0] = 1.0;
256  w[1] = 1.0;
257  }
258  else
259  {
260  int i;
261  double fac, one = 1.0, apb = alpha + beta, two = 2.0;
262 
263  z[0] = -one;
264  z[np - 1] = one;
265  jacobz(np - 2, z + 1, alpha + one, beta + one);
266  jacobfd(np, z, w, NULL, np - 1, alpha, beta);
267 
268  fac = pow(two, apb + 1) * gammaFracGammaF(np, alpha, np, 0.0) *
269  gammaFracGammaF(np, beta, np + 1, apb);
270  fac /= (np - 1);
271 
272  for (i = 0; i < np; ++i)
273  w[i] = fac / (w[i] * w[i]);
274  w[0] *= (beta + one);
275  w[np - 1] *= (alpha + one);
276  }
277 
278  return;
279 }

References Nektar::LibUtilities::beta, gammaFracGammaF(), jacobfd(), and jacobz.

Referenced by main(), Nektar::LibUtilities::TimeIntegrationSchemeSDC::TimeIntegrationSchemeSDC(), and Nektar::LibUtilities::GaussPoints::v_CalculatePoints().

◆ zwgrjm()

void Polylib::zwgrjm ( double *  z,
double *  w,
const int  np,
const double  alpha,
const double  beta 
)

Gauss-Radau-Jacobi zeros and weights with end point at z=-1.

  • Generate np Gauss-Radau-Jacobi zeros, z, and weights,w, associated with the polynomial \((1+z) P^{\alpha,\beta+1}_{np-1}(z) \).
  • Exact for polynomials of order 2np-2 or less

Definition at line 161 of file Polylib.cpp.

163 {
164 
165  if (np == 1)
166  {
167  z[0] = 0.0;
168  w[0] = 2.0;
169  }
170  else
171  {
172  int i;
173  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
174 
175  z[0] = -one;
176  jacobz(np - 1, z + 1, alpha, beta + 1);
177  jacobfd(np, z, w, NULL, np - 1, alpha, beta);
178 
179  fac = pow(two, apb) * gammaFracGammaF(np, alpha, np, 0.0) *
180  gammaFracGammaF(np, beta, np + 1, apb);
181  fac /= (beta + np);
182 
183  for (i = 0; i < np; ++i)
184  w[i] = fac * (1 - z[i]) / (w[i] * w[i]);
185  w[0] *= (beta + one);
186  }
187 
188  return;
189 }

References Nektar::LibUtilities::beta, gammaFracGammaF(), jacobfd(), and jacobz.

Referenced by main(), and Nektar::LibUtilities::GaussPoints::v_CalculatePoints().

◆ zwgrjp()

void Polylib::zwgrjp ( double *  z,
double *  w,
const int  np,
const double  alpha,
const double  beta 
)

Gauss-Radau-Jacobi zeros and weights with end point at z=1.

  • Generate np Gauss-Radau-Jacobi zeros, z, and weights,w, associated with the polynomial \((1-z) P^{\alpha+1,\beta}_{np-1}(z) \).
  • Exact for polynomials of order 2np-2 or less

Definition at line 202 of file Polylib.cpp.

204 {
205 
206  if (np == 1)
207  {
208  z[0] = 0.0;
209  w[0] = 2.0;
210  }
211  else
212  {
213  int i;
214  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
215 
216  jacobz(np - 1, z, alpha + 1, beta);
217  z[np - 1] = one;
218  jacobfd(np, z, w, NULL, np - 1, alpha, beta);
219 
220  fac = pow(two, apb) * gammaFracGammaF(np, alpha, np, 0.0) *
221  gammaFracGammaF(np, beta, np + 1, apb);
222  fac /= (alpha + np);
223 
224  for (i = 0; i < np; ++i)
225  w[i] = fac * (1 + z[i]) / (w[i] * w[i]);
226  w[np - 1] *= (alpha + one);
227  }
228 
229  return;
230 }

References Nektar::LibUtilities::beta, gammaFracGammaF(), jacobfd(), and jacobz.

Referenced by main(), Nektar::LibUtilities::TimeIntegrationSchemeSDC::TimeIntegrationSchemeSDC(), and Nektar::LibUtilities::GaussPoints::v_CalculatePoints().

◆ zwlk()

void Polylib::zwlk ( double *  z,
double *  w,
const int  npt,
const double  alpha,
const double  beta 
)

Gauss-Lobatto-Kronrod-Jacobi zeros and weights.

  • Generate npt=2*np-1 Lobatto-Kronrod Jacobi zeros, z, and weights,w, associated with the Jacobi polynomial \( P^{\alpha,\beta}_{np}(z) \),

Definition at line 468 of file Polylib.cpp.

470 {
471 
472  int np = (npt + 1) / 2;
473 
474  if (np < 4)
475  {
476  fprintf(stderr, "too few points in formula\n");
477  return;
478  }
479 
480  double endl = -1;
481  double endr = 1;
482  int i, j;
483 
484  // number of kronrod points associated with the np gauss rule
485  int kpoints = 2 * np - 1;
486 
487  // Define the number of required recurrence coefficents
488  int ncoeffs = (int)ceil(3.0 * np / 2) - 1;
489 
490  // Define arrays for the recurrence coefficients
491  double *a = new double[ncoeffs + 1];
492  double *b = new double[ncoeffs + 1];
493 
494  // Initialize a and b to zero
495  for (i = 0; i < ncoeffs + 1; i++)
496  {
497  a[i] = 0.0;
498  b[i] = 0.0;
499  }
500 
501  // Call the routine to calculate the recurrence coefficients
502  RecCoeff(ncoeffs, a, b, alpha, beta);
503 
504  double *a0 = new double[ncoeffs];
505  double *b0 = new double[ncoeffs];
506 
507  chri1(ncoeffs, a, b, a0, b0, endl);
508 
509  double *a1 = new double[ncoeffs - 1];
510  double *b1 = new double[ncoeffs - 1];
511 
512  chri1(ncoeffs - 1, a0, b0, a1, b1, endr);
513 
514  double s = b1[0] / fabs(b1[0]);
515  b1[0] = s * b1[0];
516 
517  // Finding the 2*np-1 gauss-kronrod points
518  double *z1 = new double[2 * np - 3];
519  double *w1 = new double[2 * np - 3];
520  for (i = 0; i < ncoeffs; i++)
521  {
522  z1[i] = a1[i];
523  w1[i] = b1[i];
524  }
525  JKMatrix(np - 2, z1, w1);
526  // Set up the identity matrix
527  double **zmatrix = new double *[2 * np - 3];
528  for (i = 0; i < 2 * np - 3; i++)
529  {
530  zmatrix[i] = new double[2 * np - 3];
531  for (j = 0; j < 2 * np - 3; j++)
532  {
533  zmatrix[i][j] = 0.0;
534  }
535  }
536  for (i = 0; i < 2 * np - 3; i++)
537  {
538  zmatrix[i][i] = 1.0;
539  }
540 
541  // Calculate the points and weights
542  TriQL(2 * np - 3, z1, w1, zmatrix);
543 
544  double sumW1 = 0.0;
545  double sumW1Z1 = 0.0;
546  for (i = 0; i < 2 * np - 3; i++)
547  {
548  w1[i] = s * w1[i] / (z1[i] - endl) / (z1[i] - endr);
549  sumW1 += w1[i];
550  sumW1Z1 += z1[i] * w1[i];
551  }
552 
553  double c0 = b[0] - sumW1;
554  double c1 = a[0] * b[0] - sumW1Z1;
555 
556  z[0] = endl;
557  z[2 * np - 2] = endr;
558  w[0] = (c0 * endr - c1) / (endr - endl);
559  w[2 * np - 2] = (c1 - c0 * endl) / (endr - endl);
560 
561  for (i = 1; i < kpoints - 1; i++)
562  {
563  z[i] = z1[i - 1];
564  w[i] = w1[i - 1];
565  }
566  delete[] a;
567  delete[] b;
568  delete[] a0;
569  delete[] b0;
570  delete[] a1;
571  delete[] b1;
572  delete[] z1;
573  delete[] w1;
574  for (i = 0; i < 2 * np - 3; i++)
575  {
576  delete[] zmatrix[i];
577  }
578  delete[] zmatrix;
579 }
Polylib::chri1
void chri1(int, double *, double *, double *, double *, double)
Given a weight function through the first n+1 coefficients a and b of its orthogonal polynomials thi...
Definition: Polylib.cpp:1797

References Nektar::LibUtilities::beta, chri1(), JKMatrix(), RecCoeff(), and TriQL().

Referenced by Nektar::LibUtilities::GaussPoints::v_CalculatePoints().

◆ zwrk()

void Polylib::zwrk ( double *  z,
double *  w,
const int  npt,
const double  alpha,
const double  beta 
)

Gauss-Radau-Kronrod-Jacobi zeros and weights.

  • Generate npt=2*np Radau-Kronrod Jacobi zeros, z, and weights,w, associated with the Jacobi polynomial \( P^{\alpha,\beta}_{np}(z) \),

Definition at line 362 of file Polylib.cpp.

364 {
365 
366  int np = npt / 2;
367 
368  if (np < 2)
369  {
370  fprintf(stderr, "too few points in formula\n");
371  return;
372  }
373 
374  double end0 = -1;
375 
376  int i, j;
377 
378  // number of kronrod points associated with the np gauss rule
379  int kpoints = 2 * np;
380 
381  // Define the number of required recurrence coefficents
382  int ncoeffs = (int)ceil(3.0 * np / 2);
383 
384  // Define arrays for the recurrence coefficients
385  double *a = new double[ncoeffs + 1];
386  double *b = new double[ncoeffs + 1];
387 
388  // Initialize a and b to zero
389  for (i = 0; i < ncoeffs + 1; i++)
390  {
391  a[i] = 0.0;
392  b[i] = 0.0;
393  }
394 
395  // Call the routine to calculate the recurrence coefficients
396  RecCoeff(ncoeffs, a, b, alpha, beta);
397 
398  double *a0 = new double[ncoeffs];
399  double *b0 = new double[ncoeffs];
400 
401  chri1(ncoeffs, a, b, a0, b0, end0);
402 
403  double s = b0[0] / fabs(b0[0]);
404  b0[0] = s * b0[0];
405 
406  // Finding the 2*np-1 gauss-kronrod points
407  double *z1 = new double[2 * np - 1];
408  double *w1 = new double[2 * np - 1];
409  for (i = 0; i < ncoeffs; i++)
410  {
411  z1[i] = a0[i];
412  w1[i] = b0[i];
413  }
414  JKMatrix(np - 1, z1, w1);
415  // Set up the identity matrix
416  double **zmatrix = new double *[2 * np - 1];
417  for (i = 0; i < 2 * np - 1; i++)
418  {
419  zmatrix[i] = new double[2 * np - 1];
420  for (j = 0; j < 2 * np - 1; j++)
421  {
422  zmatrix[i][j] = 0.0;
423  }
424  }
425  for (i = 0; i < 2 * np - 1; i++)
426  {
427  zmatrix[i][i] = 1.0;
428  }
429 
430  // Calculate the points and weights
431  TriQL(2 * np - 1, z1, w1, zmatrix);
432 
433  double sumW1 = 0.0;
434  for (i = 0; i < 2 * np - 1; i++)
435  {
436  w1[i] = s * w1[i] / (z1[i] - end0);
437  sumW1 += w1[i];
438  }
439 
440  z[0] = end0;
441  w[0] = b[0] - sumW1;
442  for (i = 1; i < kpoints; i++)
443  {
444  z[i] = z1[i - 1];
445  w[i] = w1[i - 1];
446  }
447 
448  delete[] a;
449  delete[] b;
450  delete[] a0;
451  delete[] b0;
452  delete[] z1;
453  delete[] w1;
454  for (i = 0; i < 2 * np - 1; i++)
455  {
456  delete[] zmatrix[i];
457  }
458  delete[] zmatrix;
459 }

References Nektar::LibUtilities::beta, chri1(), JKMatrix(), RecCoeff(), and TriQL().

Referenced by Nektar::LibUtilities::GaussPoints::v_CalculatePoints().

Generated on Sat Mar 25 2023 12:36:01 for Nektar++ by doxygen 1.9.1