Actual source code: ex6.c
1: /*$Id: ex6.c,v 1.19 2001/08/07 21:30:37 bsmith Exp $*/
3: static char help[] = "Creates a matrix using 9 pt stensil, and uses it to test MatIncreaseOverlap (needed for aditive schwarts preconditioner. \n\
4: -m <size> : problem size\n\
5: -x1, -x2 <size> : no of subdomains in x and y directions\n\n";
7: #include petscsles.h
11: int FormElementStiffness(PetscReal H,PetscScalar *Ke)
12: {
13: Ke[0] = H/6.0; Ke[1] = -.125*H; Ke[2] = H/12.0; Ke[3] = -.125*H;
14: Ke[4] = -.125*H; Ke[5] = H/6.0; Ke[6] = -.125*H; Ke[7] = H/12.0;
15: Ke[8] = H/12.0; Ke[9] = -.125*H; Ke[10] = H/6.0; Ke[11] = -.125*H;
16: Ke[12] = -.125*H; Ke[13] = H/12.0; Ke[14] = -.125*H; Ke[15] = H/6.0;
17: return 0;
18: }
21: int FormElementRhs(PetscReal x,PetscReal y,PetscReal H,PetscScalar *r)
22: {
23: r[0] = 0.; r[1] = 0.; r[2] = 0.; r[3] = 0.0;
24: return 0;
25: }
29: int main(int argc,char **args)
30: {
31: Mat C;
32: int i,m = 2,N,M,ierr,idx[4],Nsub1,Nsub2,ol=1,x1,x2;
33: PetscScalar Ke[16];
34: PetscReal x,y,h;
35: IS *is1,*is2;
36: PetscTruth flg;
38: PetscInitialize(&argc,&args,(char *)0,help);
39: PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_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(PETSC_NULL,"-x1",&x1,PETSC_NULL);
46: PetscOptionsGetInt(PETSC_NULL,"-x2",&x2,PETSC_NULL);
47: /* create stiffness matrix */
48: MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,9,PETSC_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);
64: for (ol=0; ol<m+2; ++ol) {
66: PCASMCreateSubdomains2D(m+1,m+1,x1,x2,1,0,&Nsub1,&is1);
67: MatIncreaseOverlap(C,Nsub1,is1,ol);
68: PCASMCreateSubdomains2D(m+1,m+1,x1,x2,1,ol,&Nsub2,&is2);
69:
70: PetscPrintf(PETSC_COMM_SELF,"flg == 1 => both index sets are same\n");
71: if(Nsub1 != Nsub2){
72: PetscPrintf(PETSC_COMM_SELF,"Error: No of indes sets don't match\n");
73: }
74:
75: for (i=0; i<Nsub1; ++i) {
76: ISEqual(is1[i],is2[i],&flg);
77: PetscPrintf(PETSC_COMM_SELF,"i = %d,flg = %d \n",i,flg);
78:
79: }
80: for (i=0; i<Nsub1; ++i) {ISDestroy(is1[i]);}
81: for (i=0; i<Nsub2; ++i) {ISDestroy(is2[i]);}
82:
84: PetscFree(is1);
85: PetscFree(is2);
86: }
87: MatDestroy(C);
88: PetscFinalize();
89: return 0;
90: }