Actual source code: ex15.c
1: const char help[] = "Test PetscDTSimplexQuadrature()";
3: #include <petscdt.h>
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;
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;
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 %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));
34: 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);
35: }
36: }
37: PetscCall(PetscFree(eval));
38: PetscCall(PetscQuadratureDestroy(&quad));
39: PetscFunctionReturn(PETSC_SUCCESS);
40: }
42: int main(int argc, char **argv)
43: {
44: const PetscInt dimdeg[] = {0, 20, 20, 20};
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};
52: for (PetscInt t = 0; t < 3; t++) PetscCall(testQuadrature(dim, deg, types[t]));
53: }
54: }
55: PetscCall(PetscFinalize());
56: return 0;
57: }
59: /*TEST
61: test:
62: suffix: 0
64: TEST*/