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 PetscInt num_points; 9 const PetscReal *points; 10 const PetscReal *weights; 11 PetscInt p_degree = (degree + 1) / 2; 12 PetscInt p_degree_min = degree - p_degree; 13 PetscInt Nb, Nb_min; 14 PetscReal *eval; 15 PetscQuadrature quad; 16 17 PetscFunctionBegin; 18 PetscCall(PetscDTSimplexQuadrature(dim, degree, type, &quad)); 19 PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &num_points, &points, &weights)); 20 PetscCall(PetscDTBinomialInt(dim + p_degree, dim, &Nb)); 21 PetscCall(PetscDTBinomialInt(dim + p_degree_min, dim, &Nb_min)); 22 PetscCall(PetscMalloc1(num_points * Nb, &eval)); 23 PetscCall(PetscDTPKDEvalJet(dim, num_points, points, p_degree, 0, eval)); 24 for (PetscInt i = 0; i < Nb_min; i++) { 25 for (PetscInt j = i; j < Nb; j++) { 26 PetscReal I_exact = (i == j) ? 1.0 : 0.0; 27 PetscReal I_quad = 0.0; 28 PetscReal err; 29 30 for (PetscInt q = 0; q < num_points; q++) { I_quad += weights[q] * eval[i * num_points + q] * eval[j * num_points + q]; } 31 err = PetscAbsReal(I_exact - I_quad); 32 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)); 33 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); 34 } 35 } 36 PetscCall(PetscFree(eval)); 37 PetscCall(PetscQuadratureDestroy(&quad)); 38 PetscFunctionReturn(0); 39 } 40 41 int main(int argc, char **argv) { 42 const PetscInt dimdeg[] = {0, 20, 20, 20}; 43 44 PetscFunctionBeginUser; 45 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 46 for (PetscInt dim = 0; dim <= 3; dim++) { 47 for (PetscInt deg = 0; deg <= dimdeg[dim]; deg++) { 48 const PetscDTSimplexQuadratureType types[] = {PETSCDTSIMPLEXQUAD_DEFAULT, PETSCDTSIMPLEXQUAD_CONIC, PETSCDTSIMPLEXQUAD_MINSYM}; 49 50 for (PetscInt t = 0; t < 3; t++) { PetscCall(testQuadrature(dim, deg, types[t])); } 51 } 52 } 53 PetscCall(PetscFinalize()); 54 return 0; 55 } 56 57 /*TEST 58 59 test: 60 suffix: 0 61 62 TEST*/ 63