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
testQuadrature(PetscInt dim,PetscInt degree,PetscDTSimplexQuadratureType type)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
main(int argc,char ** argv)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