1 const char help[] = "Tests PetscDTPKDEvalJet()"; 2 3 #include <petscdt.h> 4 #include <petscblaslapack.h> 5 6 static PetscErrorCode testOrthogonality(PetscInt dim, PetscInt deg) { 7 PetscQuadrature q; 8 const PetscReal *points, *weights; 9 PetscInt Npoly, npoints, i, j, k; 10 PetscReal *p; 11 12 PetscFunctionBegin; 13 PetscCall(PetscDTStroudConicalQuadrature(dim, 1, deg + 1, -1., 1., &q)); 14 PetscCall(PetscQuadratureGetData(q, NULL, NULL, &npoints, &points, &weights)); 15 PetscCall(PetscDTBinomialInt(dim + deg, dim, &Npoly)); 16 PetscCall(PetscMalloc1(Npoly * npoints, &p)); 17 PetscCall(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.; 22 23 for (k = 0; k < npoints; k++) integral += weights[k] * p[i * npoints + k] * p[j * npoints + k]; 24 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); 25 } 26 } 27 PetscCall(PetscFree(p)); 28 PetscCall(PetscQuadratureDestroy(&q)); 29 PetscFunctionReturn(0); 30 } 31 32 static PetscErrorCode testDerivativesLegendre(PetscInt dim, PetscInt deg, PetscInt k) { 33 PetscInt Np, Nk, i, j, l, d, npoints; 34 PetscRandom rand; 35 PetscReal *point; 36 PetscReal *lgndre_coeffs; 37 PetscReal *pkd_coeffs; 38 PetscReal *proj; 39 PetscReal **B; 40 PetscQuadrature q; 41 PetscReal *points1d; 42 PetscInt *degrees; 43 PetscInt *degtup, *ktup; 44 const PetscReal *points; 45 const PetscReal *weights; 46 PetscReal *lgndre_jet; 47 PetscReal **D; 48 PetscReal *pkd_jet, *pkd_jet_basis; 49 50 PetscFunctionBegin; 51 PetscCall(PetscDTBinomialInt(dim + deg, dim, &Np)); 52 PetscCall(PetscDTBinomialInt(dim + k, dim, &Nk)); 53 54 /* create the projector (because it is an orthonormal basis, the projector is the moment integrals) */ 55 PetscCall(PetscDTStroudConicalQuadrature(dim, 1, deg + 1, -1., 1., &q)); 56 PetscCall(PetscQuadratureGetData(q, NULL, NULL, &npoints, &points, &weights)); 57 PetscCall(PetscMalloc1(npoints * Np, &proj)); 58 PetscCall(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]; 61 62 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rand)); 63 PetscCall(PetscRandomSetInterval(rand, -1., 1.)); 64 65 /* create a random coefficient vector */ 66 PetscCall(PetscMalloc2(Np, &lgndre_coeffs, Np, &pkd_coeffs)); 67 for (i = 0; i < Np; i++) { PetscCall(PetscRandomGetValueReal(rand, &lgndre_coeffs[i])); } 68 69 PetscCall(PetscMalloc2(dim, °tup, dim, &ktup)); 70 PetscCall(PetscMalloc1(deg + 1, °rees)); 71 for (i = 0; i < deg + 1; i++) degrees[i] = i; 72 73 /* project the lgndre_coeffs to pkd_coeffs */ 74 PetscCall(PetscArrayzero(pkd_coeffs, Np)); 75 PetscCall(PetscMalloc1(npoints, &points1d)); 76 PetscCall(PetscMalloc1(dim, &B)); 77 for (d = 0; d < dim; d++) { 78 PetscCall(PetscMalloc1((deg + 1) * npoints, &(B[d]))); 79 /* get this coordinate */ 80 for (i = 0; i < npoints; i++) points1d[i] = points[i * dim + d]; 81 PetscCall(PetscDTLegendreEval(npoints, points1d, deg + 1, degrees, B[d], NULL, NULL)); 82 } 83 PetscCall(PetscFree(points1d)); 84 for (i = 0; i < npoints; i++) { 85 PetscReal val = 0.; 86 87 for (j = 0; j < Np; j++) { 88 PetscReal mul = lgndre_coeffs[j]; 89 PetscReal valj = 1.; 90 91 PetscCall(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++) { PetscCall(PetscFree(B[i])); } 98 PetscCall(PetscFree(B)); 99 100 /* create a random point in the biunit simplex */ 101 PetscCall(PetscMalloc1(dim, &point)); 102 for (i = 0; i < dim; i++) { PetscCall(PetscRandomGetValueReal(rand, &point[i])); } 103 for (i = dim - 1; i > 0; i--) { 104 PetscReal val = point[i]; 105 PetscInt j; 106 107 for (j = 0; j < i; j++) { point[j] = (point[j] + 1.) * (1. - val) * 0.5 - 1.; } 108 } 109 110 PetscCall(PetscMalloc3(Nk * Np, &pkd_jet_basis, Nk, &lgndre_jet, Nk, &pkd_jet)); 111 /* evaluate the jet at the point with PKD polynomials */ 112 PetscCall(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 } 118 119 /* evaluate the 1D jets of the Legendre polynomials */ 120 PetscCall(PetscMalloc1(dim, &D)); 121 for (i = 0; i < dim; i++) { 122 PetscCall(PetscMalloc1((deg + 1) * (k + 1), &(D[i]))); 123 PetscCall(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]; 129 130 PetscCall(PetscDTIndexToGradedOrder(dim, i, degtup)); 131 for (j = 0; j < Nk; j++) { 132 PetscReal val = 1.; 133 134 PetscCall(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++) { PetscCall(PetscFree(D[i])); } 140 PetscCall(PetscFree(D)); 141 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; 146 147 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); 148 } 149 150 PetscCall(PetscFree2(degtup, ktup)); 151 PetscCall(PetscFree(degrees)); 152 PetscCall(PetscFree3(pkd_jet_basis, lgndre_jet, pkd_jet)); 153 PetscCall(PetscFree(point)); 154 PetscCall(PetscFree2(lgndre_coeffs, pkd_coeffs)); 155 PetscCall(PetscFree(proj)); 156 PetscCall(PetscRandomDestroy(&rand)); 157 PetscCall(PetscQuadratureDestroy(&q)); 158 PetscFunctionReturn(0); 159 } 160 161 int main(int argc, char **argv) { 162 PetscInt dim, deg, k; 163 164 PetscFunctionBeginUser; 165 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 166 dim = 3; 167 deg = 4; 168 k = 3; 169 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for PetscDTPKDEval() tests", "none"); 170 PetscCall(PetscOptionsInt("-dim", "Dimension of the simplex", "ex9.c", dim, &dim, NULL)); 171 PetscCall(PetscOptionsInt("-degree", "The degree of the polynomial space", "ex9.c", deg, °, NULL)); 172 PetscCall(PetscOptionsInt("-k", "The number of derivatives to use in the taylor test", "ex9.c", k, &k, NULL)); 173 PetscOptionsEnd(); 174 PetscCall(testOrthogonality(dim, deg)); 175 PetscCall(testDerivativesLegendre(dim, deg, k)); 176 PetscCall(PetscFinalize()); 177 return 0; 178 } 179 180 /*TEST 181 182 test: 183 args: -dim {{1 2 3 4}} 184 185 TEST*/ 186