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

 24:       for (k = 0; k < npoints; k++) integral += weights[k] * p[i * npoints + k] * p[j * npoints + k];
 25:       PetscCheck(PetscAbsReal(integral - exact) <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "<P[%" PetscInt_FMT "], P[%" PetscInt_FMT "]> = %g != delta_{%" PetscInt_FMT ",%" PetscInt_FMT "}", i, j, (double)integral, i, j);
 26:     }
 27:   }
 28:   PetscCall(PetscFree(p));
 29:   PetscCall(PetscQuadratureDestroy(&q));
 30:   PetscFunctionReturn(PETSC_SUCCESS);
 31: }

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

 52:   PetscFunctionBegin;
 53:   PetscCall(PetscDTBinomialInt(dim + deg, dim, &Np));
 54:   PetscCall(PetscDTBinomialInt(dim + k, dim, &Nk));

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

 64:   PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rand));
 65:   PetscCall(PetscRandomSetInterval(rand, -1., 1.));

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

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

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

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

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

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

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

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

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

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

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

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

149:     PetscCheck(PetscAbsReal(diff) <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Jet mismatch between PKD and tensor Legendre bases: error %g at tolerance %g", (double)diff, (double)tol);
150:   }

152:   PetscCall(PetscFree2(degtup, ktup));
153:   PetscCall(PetscFree(degrees));
154:   PetscCall(PetscFree3(pkd_jet_basis, lgndre_jet, pkd_jet));
155:   PetscCall(PetscFree(point));
156:   PetscCall(PetscFree2(lgndre_coeffs, pkd_coeffs));
157:   PetscCall(PetscFree(proj));
158:   PetscCall(PetscRandomDestroy(&rand));
159:   PetscCall(PetscQuadratureDestroy(&q));
160:   PetscFunctionReturn(PETSC_SUCCESS);
161: }

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

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

183: /*TEST

185:   test:
186:     args: -dim {{1 2 3 4}}

188: TEST*/