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*/