Actual source code: ex9.c

  1: const char help[] = "Tests PetscDTPKDEvalJet()";

  3: #include <petscdt.h>
  4: #include <petscblaslapack.h>

  6: static PetscErrorCode testOrthogonality(PetscInt dim, PetscInt deg)
  7: {
  8:   PetscQuadrature  q;
  9:   const PetscReal *points, *weights;
 10:   PetscInt         Npoly, npoints, i, j, k;
 11:   PetscReal       *p;

 13:   PetscDTStroudConicalQuadrature(dim, 1, deg + 1, -1., 1., &q);
 14:   PetscQuadratureGetData(q, NULL, NULL, &npoints, &points, &weights);
 15:   PetscDTBinomialInt(dim + deg, dim, &Npoly);
 16:   PetscMalloc1(Npoly * npoints, &p);
 17:   PetscDTPKDEvalJet(dim, npoints, points, deg, 0, p);
 18:   for (i = 0; i < Npoly; i++) {
 19:     for (j = i; j < Npoly; j++) {
 20:       PetscReal integral = 0.;
 21:       PetscReal exact    = (i == j) ? 1. : 0.;

 23:       for (k = 0; k < npoints; k++) integral += weights[k] * p[i * npoints + k] * p[j * npoints + k];
 25:     }
 26:   }
 27:   PetscFree(p);
 28:   PetscQuadratureDestroy(&q);
 29:   return 0;
 30: }

 32: static PetscErrorCode testDerivativesLegendre(PetscInt dim, PetscInt deg, PetscInt k)
 33: {
 34:   PetscInt         Np, Nk, i, j, l, d, npoints;
 35:   PetscRandom      rand;
 36:   PetscReal       *point;
 37:   PetscReal       *lgndre_coeffs;
 38:   PetscReal       *pkd_coeffs;
 39:   PetscReal       *proj;
 40:   PetscReal      **B;
 41:   PetscQuadrature  q;
 42:   PetscReal       *points1d;
 43:   PetscInt        *degrees;
 44:   PetscInt        *degtup, *ktup;
 45:   const PetscReal *points;
 46:   const PetscReal *weights;
 47:   PetscReal       *lgndre_jet;
 48:   PetscReal      **D;
 49:   PetscReal       *pkd_jet, *pkd_jet_basis;

 51:   PetscDTBinomialInt(dim + deg, dim, &Np);
 52:   PetscDTBinomialInt(dim + k, dim, &Nk);

 54:   /* create the projector (because it is an orthonormal basis, the projector is the moment integrals) */
 55:   PetscDTStroudConicalQuadrature(dim, 1, deg + 1, -1., 1., &q);
 56:   PetscQuadratureGetData(q, NULL, NULL, &npoints, &points, &weights);
 57:   PetscMalloc1(npoints * Np, &proj);
 58:   PetscDTPKDEvalJet(dim, npoints, points, deg, 0, proj);
 59:   for (i = 0; i < Np; i++)
 60:     for (j = 0; j < npoints; j++) proj[i * npoints + j] *= weights[j];

 62:   PetscRandomCreate(PETSC_COMM_SELF, &rand);
 63:   PetscRandomSetInterval(rand, -1., 1.);

 65:   /* create a random coefficient vector */
 66:   PetscMalloc2(Np, &lgndre_coeffs, Np, &pkd_coeffs);
 67:   for (i = 0; i < Np; i++) PetscRandomGetValueReal(rand, &lgndre_coeffs[i]);

 69:   PetscMalloc2(dim, &degtup, dim, &ktup);
 70:   PetscMalloc1(deg + 1, &degrees);
 71:   for (i = 0; i < deg + 1; i++) degrees[i] = i;

 73:   /* project the lgndre_coeffs to pkd_coeffs */
 74:   PetscArrayzero(pkd_coeffs, Np);
 75:   PetscMalloc1(npoints, &points1d);
 76:   PetscMalloc1(dim, &B);
 77:   for (d = 0; d < dim; d++) {
 78:     PetscMalloc1((deg + 1) * npoints, &(B[d]));
 79:     /* get this coordinate */
 80:     for (i = 0; i < npoints; i++) points1d[i] = points[i * dim + d];
 81:     PetscDTLegendreEval(npoints, points1d, deg + 1, degrees, B[d], NULL, NULL);
 82:   }
 83:   PetscFree(points1d);
 84:   for (i = 0; i < npoints; i++) {
 85:     PetscReal val = 0.;

 87:     for (j = 0; j < Np; j++) {
 88:       PetscReal mul  = lgndre_coeffs[j];
 89:       PetscReal valj = 1.;

 91:       PetscDTIndexToGradedOrder(dim, j, degtup);
 92:       for (l = 0; l < dim; l++) valj *= B[l][i * (deg + 1) + degtup[l]];
 93:       val += mul * valj;
 94:     }
 95:     for (j = 0; j < Np; j++) pkd_coeffs[j] += proj[j * npoints + i] * val;
 96:   }
 97:   for (i = 0; i < dim; i++) PetscFree(B[i]);
 98:   PetscFree(B);

100:   /* create a random point in the biunit simplex */
101:   PetscMalloc1(dim, &point);
102:   for (i = 0; i < dim; i++) PetscRandomGetValueReal(rand, &point[i]);
103:   for (i = dim - 1; i > 0; i--) {
104:     PetscReal val = point[i];
105:     PetscInt  j;

107:     for (j = 0; j < i; j++) point[j] = (point[j] + 1.) * (1. - val) * 0.5 - 1.;
108:   }

110:   PetscMalloc3(Nk * Np, &pkd_jet_basis, Nk, &lgndre_jet, Nk, &pkd_jet);
111:   /* evaluate the jet at the point with PKD polynomials */
112:   PetscDTPKDEvalJet(dim, 1, point, deg, k, pkd_jet_basis);
113:   for (i = 0; i < Nk; i++) {
114:     PetscReal val = 0.;
115:     for (j = 0; j < Np; j++) val += pkd_coeffs[j] * pkd_jet_basis[j * Nk + i];
116:     pkd_jet[i] = val;
117:   }

119:   /* evaluate the 1D jets of the Legendre polynomials */
120:   PetscMalloc1(dim, &D);
121:   for (i = 0; i < dim; i++) {
122:     PetscMalloc1((deg + 1) * (k + 1), &(D[i]));
123:     PetscDTJacobiEvalJet(0., 0., 1, &(point[i]), deg, k, D[i]);
124:   }
125:   /* compile the 1D Legendre jets into the tensor Legendre jet */
126:   for (j = 0; j < Nk; j++) lgndre_jet[j] = 0.;
127:   for (i = 0; i < Np; i++) {
128:     PetscReal mul = lgndre_coeffs[i];

130:     PetscDTIndexToGradedOrder(dim, i, degtup);
131:     for (j = 0; j < Nk; j++) {
132:       PetscReal val = 1.;

134:       PetscDTIndexToGradedOrder(dim, j, ktup);
135:       for (l = 0; l < dim; l++) val *= D[l][degtup[l] * (k + 1) + ktup[l]];
136:       lgndre_jet[j] += mul * val;
137:     }
138:   }
139:   for (i = 0; i < dim; i++) PetscFree(D[i]);
140:   PetscFree(D);

142:   for (i = 0; i < Nk; i++) {
143:     PetscReal diff  = lgndre_jet[i] - pkd_jet[i];
144:     PetscReal scale = 1. + PetscAbsReal(lgndre_jet[i]) + PetscAbsReal(pkd_jet[i]);
145:     PetscReal tol   = 10. * PETSC_SMALL * scale;

148:   }

150:   PetscFree2(degtup, ktup);
151:   PetscFree(degrees);
152:   PetscFree3(pkd_jet_basis, lgndre_jet, pkd_jet);
153:   PetscFree(point);
154:   PetscFree2(lgndre_coeffs, pkd_coeffs);
155:   PetscFree(proj);
156:   PetscRandomDestroy(&rand);
157:   PetscQuadratureDestroy(&q);
158:   return 0;
159: }

161: int main(int argc, char **argv)
162: {
163:   PetscInt dim, deg, k;

166:   PetscInitialize(&argc, &argv, NULL, help);
167:   dim = 3;
168:   deg = 4;
169:   k   = 3;
170:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for PetscDTPKDEval() tests", "none");
171:   PetscOptionsInt("-dim", "Dimension of the simplex", "ex9.c", dim, &dim, NULL);
172:   PetscOptionsInt("-degree", "The degree of the polynomial space", "ex9.c", deg, &deg, NULL);
173:   PetscOptionsInt("-k", "The number of derivatives to use in the taylor test", "ex9.c", k, &k, NULL);
174:   PetscOptionsEnd();
175:   testOrthogonality(dim, deg);
176:   testDerivativesLegendre(dim, deg, k);
177:   PetscFinalize();
178:   return 0;
179: }

181: /*TEST

183:   test:
184:     args: -dim {{1 2 3 4}}

186: TEST*/