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