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, °tup, dim, &ktup));
72: PetscCall(PetscMalloc1(deg + 1, °rees));
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, °, 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*/