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++) I_quad += weights[q] * eval[i * num_points + q] * eval[j * num_points + q]; 32 err = PetscAbsReal(I_exact - I_quad); 33 PetscCall(PetscInfo(quad, "Dimension %" PetscInt_FMT ", degree %" PetscInt_FMT ", method %s, error in <P_PKD(%" PetscInt_FMT "),P_PKD%" PetscInt_FMT "d)> = %g\n", dim, degree, PetscDTSimplexQuadratureTypes[type], i, j, (double)err)); 34 PetscCheck(err < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Dimension %" PetscInt_FMT ", degree %" PetscInt_FMT ", method %s, error in <P_PKD(%" PetscInt_FMT "),P_PKD(%" PetscInt_FMT ")> = %g", dim, degree, PetscDTSimplexQuadratureTypes[type], i, j, (double)err); 35 } 36 } 37 PetscCall(PetscFree(eval)); 38 PetscCall(PetscQuadratureDestroy(&quad)); 39 PetscFunctionReturn(PETSC_SUCCESS); 40 } 41 42 int main(int argc, char **argv) 43 { 44 const PetscInt dimdeg[] = {0, 20, 20, 20}; 45 46 PetscFunctionBeginUser; 47 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 48 for (PetscInt dim = 0; dim <= 3; dim++) { 49 for (PetscInt deg = 0; deg <= dimdeg[dim]; deg++) { 50 const PetscDTSimplexQuadratureType types[] = {PETSCDTSIMPLEXQUAD_DEFAULT, PETSCDTSIMPLEXQUAD_CONIC, PETSCDTSIMPLEXQUAD_MINSYM}; 51 52 for (PetscInt t = 0; t < 3; t++) PetscCall(testQuadrature(dim, deg, types[t])); 53 } 54 } 55 PetscCall(PetscFinalize()); 56 return 0; 57 } 58 59 /*TEST 60 61 test: 62 suffix: 0 63 output_file: output/empty.out 64 65 TEST*/ 66