1: static char help[] = "Tests 1D cell-based discretization tools.\n\n";

  3: #include <petscdt.h>

  5: int main(int argc, char **argv)
  6: {
  7:   PetscInt  i, j, degrees[1000], ndegrees, nsrc_points, ntarget_points;
  8:   PetscReal src_points[1000], target_points[1000], *R;
  9:   PetscBool flg;

 11:   PetscFunctionBeginUser;
 12:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 13:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Discretization tools test options", NULL);
 14:   {
 15:     ndegrees   = 1000;
 16:     degrees[0] = 1;
 17:     degrees[1] = 2;
 18:     degrees[2] = 3;
 19:     PetscCall(PetscOptionsIntArray("-degrees", "list of max degrees to evaluate", "", degrees, &ndegrees, &flg));
 20:     if (!flg) ndegrees = 3;

 22:     nsrc_points   = 1000;
 23:     src_points[0] = -1.;
 24:     src_points[1] = 0.;
 25:     src_points[2] = 1.;
 26:     PetscCall(PetscOptionsRealArray("-src_points", "list of points defining intervals on which to integrate", "", src_points, &nsrc_points, &flg));
 27:     if (!flg) nsrc_points = 3;

 29:     ntarget_points   = 1000;
 30:     target_points[0] = -1.;
 31:     target_points[1] = 0.;
 32:     target_points[2] = 1.;
 33:     PetscCall(PetscOptionsRealArray("-target_points", "list of points defining intervals on which to integrate", "", target_points, &ntarget_points, &flg));
 34:     if (!flg) ntarget_points = 3;
 35:   }
 36:   PetscOptionsEnd();

 38:   PetscCall(PetscMalloc1((nsrc_points - 1) * (ntarget_points - 1), &R));
 39:   for (i = 0; i < ndegrees; i++) {
 40:     PetscCall(PetscDTReconstructPoly(degrees[i], nsrc_points - 1, src_points, ntarget_points - 1, target_points, R));
 41:     for (j = 0; j < (ntarget_points - 1) * (nsrc_points - 1); j++) { /* Truncate to zero for nicer output */
 42:       if (PetscAbs(R[j]) < 10 * PETSC_MACHINE_EPSILON) R[j] = 0;
 43:     }
 44:     for (j = 0; j < ntarget_points - 1; j++) {
 45:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Degree %" PetscInt_FMT " target interval (%g,%g)\n", degrees[i], (double)target_points[j], (double)target_points[j + 1]));
 46:       PetscCall(PetscRealView(nsrc_points - 1, R + j * (nsrc_points - 1), PETSC_VIEWER_STDOUT_WORLD));
 47:     }
 48:   }
 49:   PetscCall(PetscFree(R));
 50:   PetscCall(PetscFinalize());
 51:   return 0;
 52: }

 54: /*TEST
 55:   test:
 56:     suffix: 1
 57:     args: -degrees 1,2,3 -target_points -0.3,0,.2 -src_points -1,-.3,0,.2,1
 58: TEST*/