2: static char help[] = "Creates a matrix using 9 pt stencil, and uses it to test MatIncreaseOverlap (needed for additive Schwarz preconditioner). \n\
  3:   -m <size>       : problem size\n\
  4:   -x1, -x2 <size> : no of subdomains in x and y directions\n\n";

  6: #include <petscksp.h>

  8: PetscErrorCode FormElementStiffness(PetscReal H,PetscScalar *Ke)
  9: {
 10:   Ke[0]  = H/6.0;    Ke[1]  = -.125*H; Ke[2]  = H/12.0;   Ke[3]  = -.125*H;
 11:   Ke[4]  = -.125*H;  Ke[5]  = H/6.0;   Ke[6]  = -.125*H;  Ke[7]  = H/12.0;
 12:   Ke[8]  = H/12.0;   Ke[9]  = -.125*H; Ke[10] = H/6.0;    Ke[11] = -.125*H;
 13:   Ke[12] = -.125*H;  Ke[13] = H/12.0;  Ke[14] = -.125*H;  Ke[15] = H/6.0;
 14:   return 0;
 15: }
 16: PetscErrorCode FormElementRhs(PetscReal x,PetscReal y,PetscReal H,PetscScalar *r)
 17: {
 18:   r[0] = 0.; r[1] = 0.; r[2] = 0.; r[3] = 0.0;
 19:   return 0;
 20: }

 22: int main(int argc,char **args)
 23: {
 24:   Mat            C;
 26:   PetscInt       i,m = 2,N,M,idx[4],Nsub1,Nsub2,ol=1,x1,x2;
 27:   PetscScalar    Ke[16];
 28:   PetscReal      h;
 29:   IS             *is1,*is2,*islocal1,*islocal2;
 30:   PetscBool      flg;

 32:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 33:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 34:   N    = (m+1)*(m+1); /* dimension of matrix */
 35:   M    = m*m; /* number of elements */
 36:   h    = 1.0/m;    /* mesh width */
 37:   x1   = (m+1)/2;
 38:   x2   = x1;
 39:   PetscOptionsGetInt(NULL,NULL,"-x1",&x1,NULL);
 40:   PetscOptionsGetInt(NULL,NULL,"-x2",&x2,NULL);
 41:   /* create stiffness matrix */
 42:   MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,9,NULL,&C);

 44:   /* forms the element stiffness for the Laplacian */
 45:   FormElementStiffness(h*h,Ke);
 46:   for (i=0; i<M; i++) {
 47:     /* node numbers for the four corners of element */
 48:     idx[0] = (m+1)*(i/m) + (i % m);
 49:     idx[1] = idx[0]+1; idx[2] = idx[1] + m + 1; idx[3] = idx[2] - 1;
 50:     MatSetValues(C,4,idx,4,idx,Ke,ADD_VALUES);
 51:   }
 52:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
 53:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

 55:   for (ol=0; ol<m+2; ++ol) {

 57:     PCASMCreateSubdomains2D(m+1,m+1,x1,x2,1,0,&Nsub1,&is1,&islocal1);
 58:     MatIncreaseOverlap(C,Nsub1,is1,ol);
 59:     PCASMCreateSubdomains2D(m+1,m+1,x1,x2,1,ol,&Nsub2,&is2,&islocal2);

 61:     PetscPrintf(PETSC_COMM_SELF,"flg == 1 => both index sets are same\n");
 62:     if (Nsub1 != Nsub2) {
 63:       PetscPrintf(PETSC_COMM_SELF,"Error: No of indes sets don't match\n");
 64:     }

 66:     for (i=0; i<Nsub1; ++i) {
 67:       ISEqual(is1[i],is2[i],&flg);
 68:       PetscPrintf(PETSC_COMM_SELF,"i =  %D,flg = %d \n",i,(int)flg);

 70:     }
 71:     for (i=0; i<Nsub1; ++i) {ISDestroy(&is1[i]);}
 72:     for (i=0; i<Nsub2; ++i) {ISDestroy(&is2[i]);}
 73:     for (i=0; i<Nsub1; ++i) {ISDestroy(&islocal1[i]);}
 74:     for (i=0; i<Nsub2; ++i) {ISDestroy(&islocal2[i]);}


 77:     PetscFree(is1);
 78:     PetscFree(is2);
 79:     PetscFree(islocal1);
 80:     PetscFree(islocal2);
 81:   }
 82:   MatDestroy(&C);
 83:   PetscFinalize();
 84:   return ierr;
 85: }



 89: /*TEST

 91:    test:
 92:       args: -m 7

 94: TEST*/