2: static char help[] = "Creates a matrix using 9 pt stencil, and uses it to test MatIncreaseOverlap (needed for aditive schwarts 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>

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

 28: int main(int argc,char **args)
 29: {
 30:   Mat            C;
 32:   PetscInt       i,m = 2,N,M,idx[4],Nsub1,Nsub2,ol=1,x1,x2;
 33:   PetscScalar    Ke[16];
 34:   PetscReal      x,y,h;
 35:   IS             *is1,*is2;
 36:   PetscBool      flg;

 38:   PetscInitialize(&argc,&args,(char*)0,help);
 39:   PetscOptionsGetInt(NULL,"-m",&m,NULL);
 40:   N    = (m+1)*(m+1); /* dimension of matrix */
 41:   M    = m*m; /* number of elements */
 42:   h    = 1.0/m;    /* mesh width */
 43:   x1   = (m+1)/2;
 44:   x2   = x1;
 45:   PetscOptionsGetInt(NULL,"-x1",&x1,NULL);
 46:   PetscOptionsGetInt(NULL,"-x2",&x2,NULL);
 47:   /* create stiffness matrix */
 48:   MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,9,NULL,&C);

 50:   /* forms the element stiffness for the Laplacian */
 51:   FormElementStiffness(h*h,Ke);
 52:   for (i=0; i<M; i++) {
 53:     /* location of lower left corner of element */
 54:     x = h*(i % m); y = h*(i/m);
 55:     /* node numbers for the four corners of element */
 56:     idx[0] = (m+1)*(i/m) + (i % m);
 57:     idx[1] = idx[0]+1; idx[2] = idx[1] + m + 1; idx[3] = idx[2] - 1;
 58:     MatSetValues(C,4,idx,4,idx,Ke,ADD_VALUES);
 59:   }
 60:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
 61:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

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

 65:     PCASMCreateSubdomains2D(m+1,m+1,x1,x2,1,0,&Nsub1,&is1);
 66:     MatIncreaseOverlap(C,Nsub1,is1,ol);
 67:     PCASMCreateSubdomains2D(m+1,m+1,x1,x2,1,ol,&Nsub2,&is2);

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

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

 78:     }
 79:     for (i=0; i<Nsub1; ++i) {ISDestroy(&&is1[i]);}
 80:     for (i=0; i<Nsub2; ++i) {ISDestroy(&&is2[i]);}


 83:     PetscFree(is1);
 84:     PetscFree(is2);
 85:   }
 86:   MatDestroy(&C);
 87:   PetscFinalize();
 88:   return 0;
 89: }