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 PetscErrorCode ierr; 13 14 PetscFunctionBegin; 15 ierr = PetscDTStroudConicalQuadrature(dim, 1, deg + 1, -1., 1., &q);CHKERRQ(ierr); 16 ierr = PetscQuadratureGetData(q, NULL, NULL, &npoints, &points, &weights);CHKERRQ(ierr); 17 ierr = PetscDTBinomialInt(dim + deg, dim, &Npoly);CHKERRQ(ierr); 18 ierr = PetscMalloc1(Npoly * npoints, &p);CHKERRQ(ierr); 19 ierr = PetscDTPKDEvalJet(dim, npoints, points, deg, 0, p);CHKERRQ(ierr); 20 for (i = 0; i < Npoly; i++) { 21 for (j = i; j < Npoly; j++) { 22 PetscReal integral = 0.; 23 PetscReal exact = (i == j) ? 1. : 0.; 24 25 for (k = 0; k < npoints; k++) integral += weights[k] * p[i * npoints + k] * p[j * npoints + k]; 26 if (PetscAbsReal(integral - exact) > PETSC_SMALL) SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_PLIB, "<P[%D], P[%D]> = %g != delta_{%D,%D}", i, j, (double) integral, i, j); 27 } 28 } 29 ierr = PetscFree(p);CHKERRQ(ierr); 30 ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr); 31 PetscFunctionReturn(0); 32 } 33 34 static PetscErrorCode testDerivativesLegendre(PetscInt dim, PetscInt deg, PetscInt k) 35 { 36 PetscInt Np, Nk, i, j, l, d, npoints; 37 PetscRandom rand; 38 PetscReal *point; 39 PetscReal *lgndre_coeffs; 40 PetscReal *pkd_coeffs; 41 PetscReal *proj; 42 PetscReal **B; 43 PetscQuadrature q; 44 PetscReal *points1d; 45 PetscInt *degrees; 46 PetscInt *degtup, *ktup; 47 const PetscReal *points; 48 const PetscReal *weights; 49 PetscReal *lgndre_jet; 50 PetscReal **D; 51 PetscReal *pkd_jet, *pkd_jet_basis; 52 PetscErrorCode ierr; 53 54 PetscFunctionBegin; 55 ierr = PetscDTBinomialInt(dim + deg, dim, &Np);CHKERRQ(ierr); 56 ierr = PetscDTBinomialInt(dim + k, dim, &Nk);CHKERRQ(ierr); 57 58 /* create the projector (because it is an orthonormal basis, the projector is the moment integrals) */ 59 ierr = PetscDTStroudConicalQuadrature(dim, 1, deg + 1, -1., 1., &q);CHKERRQ(ierr); 60 ierr = PetscQuadratureGetData(q, NULL, NULL, &npoints, &points, &weights);CHKERRQ(ierr); 61 ierr = PetscMalloc1(npoints * Np, &proj);CHKERRQ(ierr); 62 ierr = PetscDTPKDEvalJet(dim, npoints, points, deg, 0, proj);CHKERRQ(ierr); 63 for (i = 0; i < Np; i++) for (j = 0; j < npoints; j++) proj[i * npoints + j] *= weights[j]; 64 65 ierr = PetscRandomCreate(PETSC_COMM_SELF, &rand);CHKERRQ(ierr); 66 ierr = PetscRandomSetInterval(rand, -1., 1.);CHKERRQ(ierr); 67 68 /* create a random coefficient vector */ 69 ierr = PetscMalloc2(Np, &lgndre_coeffs, Np, &pkd_coeffs);CHKERRQ(ierr); 70 for (i = 0; i < Np; i++) { 71 ierr = PetscRandomGetValueReal(rand, &lgndre_coeffs[i]);CHKERRQ(ierr); 72 } 73 74 ierr = PetscMalloc2(dim, °tup, dim, &ktup);CHKERRQ(ierr); 75 ierr = PetscMalloc1(deg + 1, °rees);CHKERRQ(ierr); 76 for (i = 0; i < deg + 1; i++) degrees[i] = i; 77 78 /* project the lgndre_coeffs to pkd_coeffs */ 79 ierr = PetscArrayzero(pkd_coeffs, Np);CHKERRQ(ierr); 80 ierr = PetscMalloc1(npoints, &points1d);CHKERRQ(ierr); 81 ierr = PetscMalloc1(dim, &B);CHKERRQ(ierr); 82 for (d = 0; d < dim; d++) { 83 ierr = PetscMalloc1((deg + 1)*npoints, &(B[d]));CHKERRQ(ierr); 84 /* get this coordinate */ 85 for (i = 0; i < npoints; i++) points1d[i] = points[i * dim + d]; 86 ierr = PetscDTLegendreEval(npoints, points1d, deg + 1, degrees, B[d], NULL, NULL);CHKERRQ(ierr); 87 } 88 ierr = PetscFree(points1d);CHKERRQ(ierr); 89 for (i = 0; i < npoints; i++) { 90 PetscReal val = 0.; 91 92 for (j = 0; j < Np; j++) { 93 PetscReal mul = lgndre_coeffs[j]; 94 PetscReal valj = 1.; 95 96 ierr = PetscDTIndexToGradedOrder(dim, j, degtup);CHKERRQ(ierr); 97 for (l = 0; l < dim; l++) { 98 valj *= B[l][i * (deg + 1) + degtup[l]]; 99 } 100 val += mul * valj; 101 } 102 for (j = 0; j < Np; j++) { 103 pkd_coeffs[j] += proj[j * npoints + i] * val; 104 } 105 } 106 for (i = 0; i < dim; i++) { 107 ierr = PetscFree(B[i]);CHKERRQ(ierr); 108 } 109 ierr = PetscFree(B);CHKERRQ(ierr); 110 111 /* create a random point in the biunit simplex */ 112 ierr = PetscMalloc1(dim, &point);CHKERRQ(ierr); 113 for (i = 0; i < dim; i++) { 114 ierr = PetscRandomGetValueReal(rand, &point[i]);CHKERRQ(ierr); 115 } 116 for (i = dim - 1; i > 0; i--) { 117 PetscReal val = point[i]; 118 PetscInt j; 119 120 for (j = 0; j < i; j++) { 121 point[j] = (point[j] + 1.)*(1. - val)*0.5 - 1.; 122 } 123 } 124 125 ierr = PetscMalloc3(Nk*Np, &pkd_jet_basis, Nk, &lgndre_jet, Nk, &pkd_jet);CHKERRQ(ierr); 126 /* evaluate the jet at the point with PKD polynomials */ 127 ierr = PetscDTPKDEvalJet(dim, 1, point, deg, k, pkd_jet_basis);CHKERRQ(ierr); 128 for (i = 0; i < Nk; i++) { 129 PetscReal val = 0.; 130 for (j = 0; j < Np; j++) { 131 val += pkd_coeffs[j] * pkd_jet_basis[j * Nk + i]; 132 } 133 pkd_jet[i] = val; 134 } 135 136 /* evaluate the 1D jets of the Legendre polynomials */ 137 ierr = PetscMalloc1(dim, &D);CHKERRQ(ierr); 138 for (i = 0; i < dim; i++) { 139 ierr = PetscMalloc1((deg + 1) * (k+1), &(D[i]));CHKERRQ(ierr); 140 ierr = PetscDTJacobiEvalJet(0., 0., 1, &(point[i]), deg, k, D[i]);CHKERRQ(ierr); 141 } 142 /* compile the 1D Legendre jets into the tensor Legendre jet */ 143 for (j = 0; j < Nk; j++) lgndre_jet[j] = 0.; 144 for (i = 0; i < Np; i++) { 145 PetscReal mul = lgndre_coeffs[i]; 146 147 ierr = PetscDTIndexToGradedOrder(dim, i, degtup);CHKERRQ(ierr); 148 for (j = 0; j < Nk; j++) { 149 PetscReal val = 1.; 150 151 ierr = PetscDTIndexToGradedOrder(dim, j, ktup);CHKERRQ(ierr); 152 for (l = 0; l < dim; l++) { 153 val *= D[l][degtup[l]*(k+1) + ktup[l]]; 154 } 155 lgndre_jet[j] += mul * val; 156 } 157 } 158 for (i = 0; i < dim; i++) { 159 ierr = PetscFree(D[i]);CHKERRQ(ierr); 160 } 161 ierr = PetscFree(D);CHKERRQ(ierr); 162 163 for (i = 0; i < Nk; i++) { 164 PetscReal diff = lgndre_jet[i] - pkd_jet[i]; 165 PetscReal scale = 1. + PetscAbsReal(lgndre_jet[i]) + PetscAbsReal(pkd_jet[i]); 166 PetscReal tol = 10. * PETSC_SMALL * scale; 167 168 if (PetscAbsReal(diff) > tol) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Jet mismatch between PKD and tensor Legendre bases: error %g at tolerance %g\n", (double) diff, (double) tol); 169 } 170 171 ierr = PetscFree2(degtup,ktup);CHKERRQ(ierr); 172 ierr = PetscFree(degrees);CHKERRQ(ierr); 173 ierr = PetscFree3(pkd_jet_basis, lgndre_jet, pkd_jet);CHKERRQ(ierr); 174 ierr = PetscFree(point);CHKERRQ(ierr); 175 ierr = PetscFree2(lgndre_coeffs, pkd_coeffs);CHKERRQ(ierr); 176 ierr = PetscFree(proj);CHKERRQ(ierr); 177 ierr = PetscRandomDestroy(&rand);CHKERRQ(ierr); 178 ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr); 179 PetscFunctionReturn(0); 180 } 181 182 int main(int argc, char **argv) 183 { 184 PetscInt dim, deg, k; 185 PetscErrorCode ierr; 186 187 ierr = PetscInitialize(&argc, &argv, NULL, help); if (ierr) return ierr; 188 dim = 3; 189 deg = 4; 190 k = 3; 191 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Options for PetscDTPKDEval() tests","none");CHKERRQ(ierr); 192 ierr = PetscOptionsInt("-dim", "Dimension of the simplex","ex9.c",dim,&dim,NULL);CHKERRQ(ierr); 193 ierr = PetscOptionsInt("-degree", "The degree of the polynomial space","ex9.c",deg,°,NULL);CHKERRQ(ierr); 194 ierr = PetscOptionsInt("-k", "The number of derivatives to use in the taylor test","ex9.c",k,&k,NULL);CHKERRQ(ierr); 195 ierr = PetscOptionsEnd();CHKERRQ(ierr); 196 ierr = testOrthogonality(dim, deg);CHKERRQ(ierr); 197 ierr = testDerivativesLegendre(dim, deg, k);CHKERRQ(ierr); 198 ierr = PetscFinalize(); 199 return ierr; 200 } 201 202 /*TEST 203 204 test: 205 args: -dim {{1 2 3 4}} 206 207 TEST*/ 208