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