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:   PetscDTSimplexQuadrature(dim, degree, type, &quad);
 19:   PetscQuadratureGetData(quad, NULL, NULL, &num_points, &points, &weights);
 20:   PetscDTBinomialInt(dim + p_degree, dim, &Nb);
 21:   PetscDTBinomialInt(dim + p_degree_min, dim, &Nb_min);
 22:   PetscMalloc1(num_points * Nb, &eval);
 23:   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;

 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:       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:     }
 35:   }
 36:   PetscFree(eval);
 37:   PetscQuadratureDestroy(&quad);
 38:   return 0;
 39: }

 41: int main(int argc, char **argv)
 42: {
 43:   const PetscInt dimdeg[] = {0, 20, 20, 20};

 46:   PetscInitialize(&argc, &argv, NULL, help);
 47:   for (PetscInt dim = 0; dim <= 3; dim++) {
 48:     for (PetscInt deg = 0; deg <= dimdeg[dim]; deg++) {
 49:       const PetscDTSimplexQuadratureType types[] = {PETSCDTSIMPLEXQUAD_DEFAULT, PETSCDTSIMPLEXQUAD_CONIC, PETSCDTSIMPLEXQUAD_MINSYM};

 51:       for (PetscInt t = 0; t < 3; t++) testQuadrature(dim, deg, types[t]);
 52:     }
 53:   }
 54:   PetscFinalize();
 55:   return 0;
 56: }

 58: /*TEST

 60:   test:
 61:     suffix: 0

 63: TEST*/