Actual source code: ex1.c
petsc-3.4.2 2013-07-02

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

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

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

     15:   PetscMalloc3(npoints*ndegrees,PetscReal,&B,npoints*ndegrees,PetscReal,&D,npoints*ndegrees,PetscReal,&D2);
     16:   PetscDTLegendreEval(npoints,points,ndegrees,degrees,B,D,D2);
     17:   PetscPrintf(PETSC_COMM_WORLD,"%s\n",name);
     18:   for (i=0; i<npoints; i++) {
     19:     for (j=0; j<ndegrees; j++) {
     20:       PetscPrintf(PETSC_COMM_WORLD,"degree %D at %12.4G: B=%12.4G  D=%12.4G  D2=%12.4G\n",degrees[j],points[i],B[i*ndegrees+j],D[i*ndegrees+j],D2[i*ndegrees+j]);
     21:     }
     22:   }
     23:   PetscFree3(B,D,D2);
     24:   return(0);
     25: }

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

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

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

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

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