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

  3: #include <petscdt.h>
  4: #include <petscviewer.h>

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

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

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

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

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