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 PetscCheckFalse(PetscAbsReal(integral - exact) > PETSC_SMALL,PETSC_COMM_SELF, PETSC_ERR_PLIB, "<P[%D], P[%D]> = %g != delta_{%D,%D}", i, j, (double) integral, i, j); 26 } 27 } 28 PetscCall(PetscFree(p)); 29 PetscCall(PetscQuadratureDestroy(&q)); 30 PetscFunctionReturn(0); 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++) for (j = 0; j < npoints; j++) proj[i * npoints + j] *= weights[j]; 62 63 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rand)); 64 PetscCall(PetscRandomSetInterval(rand, -1., 1.)); 65 66 /* create a random coefficient vector */ 67 PetscCall(PetscMalloc2(Np, &lgndre_coeffs, Np, &pkd_coeffs)); 68 for (i = 0; i < Np; i++) { 69 PetscCall(PetscRandomGetValueReal(rand, &lgndre_coeffs[i])); 70 } 71 72 PetscCall(PetscMalloc2(dim, °tup, dim, &ktup)); 73 PetscCall(PetscMalloc1(deg + 1, °rees)); 74 for (i = 0; i < deg + 1; i++) degrees[i] = i; 75 76 /* project the lgndre_coeffs to pkd_coeffs */ 77 PetscCall(PetscArrayzero(pkd_coeffs, Np)); 78 PetscCall(PetscMalloc1(npoints, &points1d)); 79 PetscCall(PetscMalloc1(dim, &B)); 80 for (d = 0; d < dim; d++) { 81 PetscCall(PetscMalloc1((deg + 1)*npoints, &(B[d]))); 82 /* get this coordinate */ 83 for (i = 0; i < npoints; i++) points1d[i] = points[i * dim + d]; 84 PetscCall(PetscDTLegendreEval(npoints, points1d, deg + 1, degrees, B[d], NULL, NULL)); 85 } 86 PetscCall(PetscFree(points1d)); 87 for (i = 0; i < npoints; i++) { 88 PetscReal val = 0.; 89 90 for (j = 0; j < Np; j++) { 91 PetscReal mul = lgndre_coeffs[j]; 92 PetscReal valj = 1.; 93 94 PetscCall(PetscDTIndexToGradedOrder(dim, j, degtup)); 95 for (l = 0; l < dim; l++) { 96 valj *= B[l][i * (deg + 1) + degtup[l]]; 97 } 98 val += mul * valj; 99 } 100 for (j = 0; j < Np; j++) { 101 pkd_coeffs[j] += proj[j * npoints + i] * val; 102 } 103 } 104 for (i = 0; i < dim; i++) { 105 PetscCall(PetscFree(B[i])); 106 } 107 PetscCall(PetscFree(B)); 108 109 /* create a random point in the biunit simplex */ 110 PetscCall(PetscMalloc1(dim, &point)); 111 for (i = 0; i < dim; i++) { 112 PetscCall(PetscRandomGetValueReal(rand, &point[i])); 113 } 114 for (i = dim - 1; i > 0; i--) { 115 PetscReal val = point[i]; 116 PetscInt j; 117 118 for (j = 0; j < i; j++) { 119 point[j] = (point[j] + 1.)*(1. - val)*0.5 - 1.; 120 } 121 } 122 123 PetscCall(PetscMalloc3(Nk*Np, &pkd_jet_basis, Nk, &lgndre_jet, Nk, &pkd_jet)); 124 /* evaluate the jet at the point with PKD polynomials */ 125 PetscCall(PetscDTPKDEvalJet(dim, 1, point, deg, k, pkd_jet_basis)); 126 for (i = 0; i < Nk; i++) { 127 PetscReal val = 0.; 128 for (j = 0; j < Np; j++) { 129 val += pkd_coeffs[j] * pkd_jet_basis[j * Nk + i]; 130 } 131 pkd_jet[i] = val; 132 } 133 134 /* evaluate the 1D jets of the Legendre polynomials */ 135 PetscCall(PetscMalloc1(dim, &D)); 136 for (i = 0; i < dim; i++) { 137 PetscCall(PetscMalloc1((deg + 1) * (k+1), &(D[i]))); 138 PetscCall(PetscDTJacobiEvalJet(0., 0., 1, &(point[i]), deg, k, D[i])); 139 } 140 /* compile the 1D Legendre jets into the tensor Legendre jet */ 141 for (j = 0; j < Nk; j++) lgndre_jet[j] = 0.; 142 for (i = 0; i < Np; i++) { 143 PetscReal mul = lgndre_coeffs[i]; 144 145 PetscCall(PetscDTIndexToGradedOrder(dim, i, degtup)); 146 for (j = 0; j < Nk; j++) { 147 PetscReal val = 1.; 148 149 PetscCall(PetscDTIndexToGradedOrder(dim, j, ktup)); 150 for (l = 0; l < dim; l++) { 151 val *= D[l][degtup[l]*(k+1) + ktup[l]]; 152 } 153 lgndre_jet[j] += mul * val; 154 } 155 } 156 for (i = 0; i < dim; i++) { 157 PetscCall(PetscFree(D[i])); 158 } 159 PetscCall(PetscFree(D)); 160 161 for (i = 0; i < Nk; i++) { 162 PetscReal diff = lgndre_jet[i] - pkd_jet[i]; 163 PetscReal scale = 1. + PetscAbsReal(lgndre_jet[i]) + PetscAbsReal(pkd_jet[i]); 164 PetscReal tol = 10. * PETSC_SMALL * scale; 165 166 PetscCheckFalse(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); 167 } 168 169 PetscCall(PetscFree2(degtup,ktup)); 170 PetscCall(PetscFree(degrees)); 171 PetscCall(PetscFree3(pkd_jet_basis, lgndre_jet, pkd_jet)); 172 PetscCall(PetscFree(point)); 173 PetscCall(PetscFree2(lgndre_coeffs, pkd_coeffs)); 174 PetscCall(PetscFree(proj)); 175 PetscCall(PetscRandomDestroy(&rand)); 176 PetscCall(PetscQuadratureDestroy(&q)); 177 PetscFunctionReturn(0); 178 } 179 180 int main(int argc, char **argv) 181 { 182 PetscInt dim, deg, k; 183 PetscErrorCode ierr; 184 185 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 186 dim = 3; 187 deg = 4; 188 k = 3; 189 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Options for PetscDTPKDEval() tests","none");PetscCall(ierr); 190 PetscCall(PetscOptionsInt("-dim", "Dimension of the simplex","ex9.c",dim,&dim,NULL)); 191 PetscCall(PetscOptionsInt("-degree", "The degree of the polynomial space","ex9.c",deg,°,NULL)); 192 PetscCall(PetscOptionsInt("-k", "The number of derivatives to use in the taylor test","ex9.c",k,&k,NULL)); 193 ierr = PetscOptionsEnd();PetscCall(ierr); 194 PetscCall(testOrthogonality(dim, deg)); 195 PetscCall(testDerivativesLegendre(dim, deg, k)); 196 PetscCall(PetscFinalize()); 197 return 0; 198 } 199 200 /*TEST 201 202 test: 203 args: -dim {{1 2 3 4}} 204 205 TEST*/ 206