Actual source code: ex2.c
1: static char help[] = "Tests 1D cell-based discretization tools.\n\n";
3: #include <petscdt.h>
4: #include <petscviewer.h>
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;
12: PetscFunctionBeginUser;
13: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
14: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Discretization tools test options", NULL);
15: {
16: ndegrees = 1000;
17: degrees[0] = 1;
18: degrees[1] = 2;
19: degrees[2] = 3;
20: PetscCall(PetscOptionsIntArray("-degrees", "list of max degrees to evaluate", "", degrees, &ndegrees, &flg));
21: if (!flg) ndegrees = 3;
23: nsrc_points = 1000;
24: src_points[0] = -1.;
25: src_points[1] = 0.;
26: src_points[2] = 1.;
27: PetscCall(PetscOptionsRealArray("-src_points", "list of points defining intervals on which to integrate", "", src_points, &nsrc_points, &flg));
28: if (!flg) nsrc_points = 3;
30: ntarget_points = 1000;
31: target_points[0] = -1.;
32: target_points[1] = 0.;
33: target_points[2] = 1.;
34: PetscCall(PetscOptionsRealArray("-target_points", "list of points defining intervals on which to integrate", "", target_points, &ntarget_points, &flg));
35: if (!flg) ntarget_points = 3;
36: }
37: PetscOptionsEnd();
39: PetscCall(PetscMalloc1((nsrc_points - 1) * (ntarget_points - 1), &R));
40: for (i = 0; i < ndegrees; i++) {
41: PetscCall(PetscDTReconstructPoly(degrees[i], nsrc_points - 1, src_points, ntarget_points - 1, target_points, R));
42: for (j = 0; j < (ntarget_points - 1) * (nsrc_points - 1); j++) { /* Truncate to zero for nicer output */
43: if (PetscAbs(R[j]) < 10 * PETSC_MACHINE_EPSILON) R[j] = 0;
44: }
45: for (j = 0; j < ntarget_points - 1; j++) {
46: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Degree %" PetscInt_FMT " target interval (%g,%g)\n", degrees[i], (double)target_points[j], (double)target_points[j + 1]));
47: PetscCall(PetscRealView(nsrc_points - 1, R + j * (nsrc_points - 1), PETSC_VIEWER_STDOUT_WORLD));
48: }
49: }
50: PetscCall(PetscFree(R));
51: PetscCall(PetscFinalize());
52: return 0;
53: }
55: /*TEST
56: test:
57: suffix: 1
58: args: -degrees 1,2,3 -target_points -0.3,0,.2 -src_points -1,-.3,0,.2,1
59: TEST*/