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*/