1 const char help[] = "Tests PetscDTPKDEvalJet()";
2
3 #include <petscdt.h>
4 #include <petscblaslapack.h>
5
testOrthogonality(PetscInt dim,PetscInt deg)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;
12
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.;
23
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 }
32
testDerivativesLegendre(PetscInt dim,PetscInt deg,PetscInt k)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;
51
52 PetscFunctionBegin;
53 PetscCall(PetscDTBinomialInt(dim + deg, dim, &Np));
54 PetscCall(PetscDTBinomialInt(dim + k, dim, &Nk));
55
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];
63
64 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rand));
65 PetscCall(PetscRandomSetInterval(rand, -1., 1.));
66
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]));
70
71 PetscCall(PetscMalloc2(dim, °tup, dim, &ktup));
72 PetscCall(PetscMalloc1(deg + 1, °rees));
73 for (i = 0; i < deg + 1; i++) degrees[i] = i;
74
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.;
88
89 for (j = 0; j < Np; j++) {
90 PetscReal mul = lgndre_coeffs[j];
91 PetscReal valj = 1.;
92
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));
101
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;
108
109 for (j = 0; j < i; j++) point[j] = (point[j] + 1.) * (1. - val) * 0.5 - 1.;
110 }
111
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 }
120
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];
131
132 PetscCall(PetscDTIndexToGradedOrder(dim, i, degtup));
133 for (j = 0; j < Nk; j++) {
134 PetscReal val = 1.;
135
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));
143
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;
148
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 }
151
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 }
162
main(int argc,char ** argv)163 int main(int argc, char **argv)
164 {
165 PetscInt dim, deg, k;
166
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 }
182
183 /*TEST
184
185 test:
186 args: -dim {{1 2 3 4}}
187 output_file: output/empty.out
188
189 TEST*/
190