xref: /petsc/src/dm/dt/tests/ex2.c (revision ef46b1a67e276116c83b5d4ce8efc2932ea4fc0a)
1 static char help[] = "Tests 1D cell-based discretization tools.\n\n";
2 
3 #include <petscdt.h>
4 #include <petscviewer.h>
5 
6 int main(int argc,char **argv)
7 {
8   PetscInt       i,j,degrees[1000],ndegrees,nsrc_points,ntarget_points;
9   PetscReal      src_points[1000],target_points[1000],*R;
10   PetscBool      flg;
11 
12   PetscCall(PetscInitialize(&argc,&argv,(char*)0,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;
21 
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;
28 
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();
37 
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 %D 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 }
53 
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*/
59