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,&B,npoints*ndegrees,&D,npoints*ndegrees,&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],(double)points[i],(double)B[i*ndegrees+j],(double)D[i*ndegrees+j],(double)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:     if (PetscAbs(zeroth) < 1e-10) zeroth = 0.;
 78:     if (PetscAbs(first)  < 1e-10) first  = 0.;
 79:     if (PetscAbs(second) < 1e-10) second = 0.;
 80:     PetscPrintf(PETSC_COMM_WORLD,"Moment error: zeroth=%g, first=%g, second=%g\n",(double)(-zeroth),(double)(-first),(double)(-second));
 81:   }
 82:   CheckPoints("Gauss points",npoints,points,ndegrees,degrees);
 83:   PetscFinalize();
 84:   return 0;
 85: }