1 const char help[] = "Test PetscDTSimplexQuadrature()"; 2 3 #include <petscdt.h> 4 5 // if we trust the PKD polynomials (tested in ex9), then we can see how well the quadrature integrates 6 // the mass matrix, which should be the identity 7 static PetscErrorCode testQuadrature(PetscInt dim, PetscInt degree, PetscDTSimplexQuadratureType type) 8 { 9 PetscInt num_points; 10 const PetscReal *points; 11 const PetscReal *weights; 12 PetscInt p_degree = (degree + 1) / 2; 13 PetscInt p_degree_min = degree - p_degree; 14 PetscInt Nb, Nb_min; 15 PetscReal *eval; 16 PetscQuadrature quad; 17 18 PetscFunctionBegin; 19 PetscCall(PetscDTSimplexQuadrature(dim, degree, type, &quad)); 20 PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &num_points, &points, &weights)); 21 PetscCall(PetscDTBinomialInt(dim + p_degree, dim, &Nb)); 22 PetscCall(PetscDTBinomialInt(dim + p_degree_min, dim, &Nb_min)); 23 PetscCall(PetscMalloc1(num_points * Nb, &eval)); 24 PetscCall(PetscDTPKDEvalJet(dim, num_points, points, p_degree, 0, eval)); 25 for (PetscInt i = 0; i < Nb_min; i++) { 26 for (PetscInt j = i; j < Nb; j++) { 27 PetscReal I_exact = (i == j) ? 1.0 : 0.0; 28 PetscReal I_quad = 0.0; 29 PetscReal err; 30 31 for (PetscInt q = 0; q < num_points; q++) { 32 I_quad += weights[q] * eval[i * num_points + q] * eval[j * num_points + q]; 33 } 34 err = PetscAbsReal(I_exact - I_quad); 35 PetscCall(PetscInfo(quad, "Dimension %d, degree %d, method %s, error in <P_PKD(%d),P_PKD(%d)> = %g\n", (int) dim, (int) degree, PetscDTSimplexQuadratureTypes[type], (int) i, (int) j, (double) err)); 36 PetscCheck(err < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Dimension %d, degree %d, method %s, error in <P_PKD(%d),P_PKD(%d)> = %g", (int) dim, (int) degree, PetscDTSimplexQuadratureTypes[type], (int) i, (int) j, (double) err); 37 } 38 } 39 PetscCall(PetscFree(eval)); 40 PetscCall(PetscQuadratureDestroy(&quad)); 41 PetscFunctionReturn(0); 42 } 43 44 int main(int argc, char **argv) 45 { 46 const PetscInt dimdeg[] = {0, 20, 20, 20}; 47 48 PetscFunctionBeginUser; 49 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 50 for (PetscInt dim = 0; dim <= 3; dim++) { 51 for (PetscInt deg = 0; deg <= dimdeg[dim]; deg++) { 52 const PetscDTSimplexQuadratureType types[] = {PETSCDTSIMPLEXQUAD_DEFAULT, PETSCDTSIMPLEXQUAD_CONIC, PETSCDTSIMPLEXQUAD_MINSYM}; 53 54 for (PetscInt t = 0; t < 3; t++) { 55 PetscCall(testQuadrature(dim, deg, types[t])); 56 } 57 } 58 } 59 PetscCall(PetscFinalize()); 60 return 0; 61 } 62 63 /*TEST 64 65 test: 66 suffix: 0 67 68 TEST*/ 69