2:  #include <petscgll.h>
  3:  #include <petscdraw.h>
  4:  #include <petscviewer.h>
  5:  #include <petscksp.h>

  7: /*
  8:       Solves -Laplacian u = f,  u(-1) = u(1) = 0 with a single spectral element for n = 4 to N GLL points

 10:       Plots the L_2 norm of the error (evaluated via the GLL nodes and weights) as a function of n.

 12: */
 13: PetscErrorCode ComputeSolution(PetscGLL *gll,Vec x)
 14: {
 16:   PetscInt       i,n;
 17:   PetscScalar    *xx;
 18:   PetscReal      xd;

 21:   VecGetSize(x,&n);
 22:   VecGetArray(x,&xx);
 23:   for (i=0; i<n; i++) {
 24:     xd    = gll->nodes[i];
 25:     xx[i] = (xd*xd - 1.0)*PetscCosReal(5.*PETSC_PI*xd);
 26:   }
 27:   VecRestoreArray(x,&xx);
 28:   return(0);
 29: }

 31: /*
 32:       Evaluates \integral_{-1}^{1} f*v_i  where v_i is the ith basis polynomial via the GLL nodes and weights, since the v_i
 33:       basis function is zero at all nodes except the ith one the integral is simply the weight_i * f(node_i)
 34: */
 35: PetscErrorCode ComputeRhs(PetscGLL *gll,Vec b)
 36: {
 38:   PetscInt       i,n;
 39:   PetscScalar    *bb;
 40:   PetscReal      xd;

 43:   VecGetSize(b,&n);
 44:   VecGetArray(b,&bb);
 45:   for (i=0; i<n; i++) {
 46:     xd    = gll->nodes[i];
 47:     bb[i] = -gll->weights[i]*(-20.*PETSC_PI*xd*PetscSinReal(5.*PETSC_PI*xd) + (2. - (5.*PETSC_PI)*(5.*PETSC_PI)*(xd*xd - 1.))*PetscCosReal(5.*PETSC_PI*xd));
 48:   }
 49:   VecRestoreArray(b,&bb);
 50:   return(0);
 51: }

 53: int main(int argc,char **args)
 54: {
 56:   PetscGLL       gll;
 57:   PetscInt       N = 80,n;
 58:   PetscReal      **A;
 59:   Mat            K;
 60:   KSP            ksp;
 61:   PC             pc;
 62:   Vec            x,b;
 63:   PetscInt       rows[2];
 64:   PetscReal      norm,xc,yc;
 65:   PetscScalar    *f;
 66:   PetscDraw      draw;
 67:   PetscDrawLG    lg;
 68:   PetscDrawAxis  axis;

 70:   PetscInitialize(&argc,&args,NULL,NULL);if (ierr) return ierr;
 71:   PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);

 73:   PetscDrawCreate(PETSC_COMM_SELF,NULL,"Log(Error norm) vs Number of GLL points",0,0,500,500,&draw);
 74:   PetscDrawSetFromOptions(draw);
 75:   PetscDrawLGCreate(draw,1,&lg);
 76:   PetscDrawLGSetUseMarkers(lg,PETSC_TRUE);
 77:   PetscDrawLGGetAxis(lg,&axis);
 78:   PetscDrawAxisSetLabels(axis,NULL,"Number of GLL points","Log(Error Norm)");

 80:   for (n=4; n<N; n+=2) {
 81:     /*
 82:        gll simply contains the GLL node and weight values
 83:     */
 84:     PetscGLLCreate(n,PETSCGLL_VIA_LINEARALGEBRA,&gll);
 85:     /*
 86:        Creates the element stiffness matrix for the given gll
 87:     */
 88:     PetscGLLElementLaplacianCreate(&gll,&A);
 89:     MatCreateSeqDense(PETSC_COMM_SELF,n,n,&A[0][0],&K);
 90:     rows[0] = 0;
 91:     rows[1] = n-1;
 92:     KSPCreate(PETSC_COMM_SELF,&ksp);
 93:     MatCreateVecs(K,&x,&b);
 94:     ComputeRhs(&gll,b);
 95:     /*
 96:         Replace the first and last rows/columns of the matrix with the identity to obtain the zero Dirichlet boundary conditions
 97:     */
 98:     MatZeroRowsColumns(K,2,rows,1.0,x,b);
 99:     KSPSetOperators(ksp,K,K);
100:     KSPGetPC(ksp,&pc);
101:     PCSetType(pc,PCLU);
102:     KSPSetFromOptions(ksp);
103:     KSPSolve(ksp,b,x);

105:     /* compute the error to the continium problem */
106:     ComputeSolution(&gll,b);
107:     VecAXPY(x,-1.0,b);

109:     /* compute the L^2 norm of the error */
110:     VecGetArray(x,&f);
111:     PetscGLLIntegrate(&gll,f,&norm);
112:     VecRestoreArray(x,&f);
113:     norm = PetscSqrtReal(norm);
114:     PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF,"L^2 norm of the error %D %g\n",n,(double)norm);
115:     xc   = (PetscReal)n;
116:     yc   = PetscLog10Real(norm);
117:     PetscDrawLGAddPoint(lg,&xc,&yc);
118:     PetscDrawLGDraw(lg);

120:     VecDestroy(&b);
121:     VecDestroy(&x);
122:     KSPDestroy(&ksp);
123:     MatDestroy(&K);
124:     PetscGLLElementLaplacianDestroy(&gll,&A);
125:     PetscGLLDestroy(&gll);
126:   }
127:   PetscDrawSetPause(draw,-2);
128:   PetscDrawLGDestroy(&lg);
129:   PetscDrawDestroy(&draw);
130:   PetscFinalize();
131:   return ierr;
132: }

134: /*TEST

136:   build:
137:       requires: !complex

139:    test:

141: TEST*/