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*/