Actual source code: ex6.c


  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;
 25:   PetscInt       i,m = 2,N,M,idx[4],Nsub1,Nsub2,ol=1,x1,x2;
 26:   PetscScalar    Ke[16];
 27:   PetscReal      h;
 28:   IS             *is1,*is2,*islocal1,*islocal2;
 29:   PetscBool      flg;

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

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

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

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

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

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

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

 75:     PetscFree(is1);
 76:     PetscFree(is2);
 77:     PetscFree(islocal1);
 78:     PetscFree(islocal2);
 79:   }
 80:   MatDestroy(&C);
 81:   PetscFinalize();
 82:   return 0;
 83: }

 85: /*TEST

 87:    test:
 88:       args: -m 7

 90: TEST*/