Actual source code: ex1.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

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

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

      6: static PetscErrorCode CheckPoints(const char *name,PetscInt npoints,const PetscReal *points,PetscInt ndegrees,const PetscInt *degrees)
      7: {
      9:   PetscReal      *B,*D,*D2;
     10:   PetscInt       i,j;

     13:   PetscMalloc3(npoints*ndegrees,&B,npoints*ndegrees,&D,npoints*ndegrees,&D2);
     14:   PetscDTLegendreEval(npoints,points,ndegrees,degrees,B,D,D2);
     15:   PetscPrintf(PETSC_COMM_WORLD,"%s\n",name);
     16:   for (i=0; i<npoints; i++) {
     17:     for (j=0; j<ndegrees; j++) {
     18:       PetscReal b,d,d2;
     19:       b = B[i*ndegrees+j];
     20:       d = D[i*ndegrees+j];
     21:       d2 = D2[i*ndegrees+j];
     22:       if (PetscAbsReal(b) < PETSC_SMALL) b   = 0;
     23:       if (PetscAbsReal(d) < PETSC_SMALL) d   = 0;
     24:       if (PetscAbsReal(d2) < PETSC_SMALL) d2 = 0;
     25:       PetscPrintf(PETSC_COMM_WORLD,"degree %D at %12.4g: B=%12.4g  D=%12.4g  D2=%12.4g\n",degrees[j],(double)points[i],(double)b,(double)d,(double)d2);
     26:     }
     27:   }
     28:   PetscFree3(B,D,D2);
     29:   return(0);
     30: }

     32: int main(int argc,char **argv)
     33: {
     35:   PetscInt       degrees[1000],ndegrees,npoints,two;
     36:   PetscReal      points[1000],weights[1000],interval[2];
     37:   PetscBool      flg;

     39:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
     40:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Discretization tools test options",NULL);
     41:   {
     42:     ndegrees   = 1000;
     43:     degrees[0] = 0;
     44:     degrees[1] = 1;
     45:     degrees[2] = 2;
     46:     PetscOptionsIntArray("-degrees","list of degrees to evaluate","",degrees,&ndegrees,&flg);

     48:     if (!flg) ndegrees = 3;
     49:     npoints   = 1000;
     50:     points[0] = 0.0;
     51:     points[1] = -0.5;
     52:     points[2] = 1.0;
     53:     PetscOptionsRealArray("-points","list of points at which to evaluate","",points,&npoints,&flg);

     55:     if (!flg) npoints = 3;
     56:     two         = 2;
     57:     interval[0] = -1.;
     58:     interval[1] = 1.;
     59:     PetscOptionsRealArray("-interval","interval on which to construct quadrature","",interval,&two,NULL);
     60:   }
     61:   PetscOptionsEnd();
     62:   CheckPoints("User-provided points",npoints,points,ndegrees,degrees);

     64:   PetscDTGaussQuadrature(npoints,interval[0],interval[1],points,weights);
     65:   PetscPrintf(PETSC_COMM_WORLD,"Quadrature weights\n");
     66:   PetscRealView(npoints,weights,PETSC_VIEWER_STDOUT_WORLD);
     67:   {
     68:     PetscReal a = interval[0],b = interval[1],zeroth,first,second;
     69:     PetscInt  i;
     70:     zeroth = b - a;
     71:     first  = (b*b - a*a)/2;
     72:     second = (b*b*b - a*a*a)/3;
     73:     for (i=0; i<npoints; i++) {
     74:       zeroth -= weights[i];
     75:       first  -= weights[i] * points[i];
     76:       second -= weights[i] * PetscSqr(points[i]);
     77:     }
     78:     if (PetscAbs(zeroth) < 1e-10) zeroth = 0.;
     79:     if (PetscAbs(first)  < 1e-10) first  = 0.;
     80:     if (PetscAbs(second) < 1e-10) second = 0.;
     81:     PetscPrintf(PETSC_COMM_WORLD,"Moment error: zeroth=%g, first=%g, second=%g\n",(double)(-zeroth),(double)(-first),(double)(-second));
     82:   }
     83:   CheckPoints("Gauss points",npoints,points,ndegrees,degrees);
     84:   PetscFinalize();
     85:   return ierr;
     86: }

     88: /*TEST
     89:   test:
     90:     suffix: 1
     91:     args: -degrees 1,2,3,4,5 -points 0,.2,-.5,.8,.9,1 -interval -.5,1
     92: TEST*/