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 { 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 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 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 188 TEST*/ 189