1 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"; 5 6 #include <petscksp.h> 7 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 } 21 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; 30 31 PetscFunctionBeginUser; 32 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 33 PetscCall(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 PetscCall(PetscOptionsGetInt(NULL,NULL,"-x1",&x1,NULL)); 40 PetscCall(PetscOptionsGetInt(NULL,NULL,"-x2",&x2,NULL)); 41 /* create stiffness matrix */ 42 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,9,NULL,&C)); 43 44 /* forms the element stiffness for the Laplacian */ 45 PetscCall(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 PetscCall(MatSetValues(C,4,idx,4,idx,Ke,ADD_VALUES)); 51 } 52 PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 53 PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 54 55 for (ol=0; ol<m+2; ++ol) { 56 57 PetscCall(PCASMCreateSubdomains2D(m+1,m+1,x1,x2,1,0,&Nsub1,&is1,&islocal1)); 58 PetscCall(MatIncreaseOverlap(C,Nsub1,is1,ol)); 59 PetscCall(PCASMCreateSubdomains2D(m+1,m+1,x1,x2,1,ol,&Nsub2,&is2,&islocal2)); 60 61 PetscCall(PetscPrintf(PETSC_COMM_SELF,"flg == 1 => both index sets are same\n")); 62 if (Nsub1 != Nsub2) { 63 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: No of indes sets don't match\n")); 64 } 65 66 for (i=0; i<Nsub1; ++i) { 67 PetscCall(ISEqual(is1[i],is2[i],&flg)); 68 PetscCall(PetscPrintf(PETSC_COMM_SELF,"i = %" PetscInt_FMT ",flg = %d \n",i,(int)flg)); 69 70 } 71 for (i=0; i<Nsub1; ++i) PetscCall(ISDestroy(&is1[i])); 72 for (i=0; i<Nsub2; ++i) PetscCall(ISDestroy(&is2[i])); 73 for (i=0; i<Nsub1; ++i) PetscCall(ISDestroy(&islocal1[i])); 74 for (i=0; i<Nsub2; ++i) PetscCall(ISDestroy(&islocal2[i])); 75 76 PetscCall(PetscFree(is1)); 77 PetscCall(PetscFree(is2)); 78 PetscCall(PetscFree(islocal1)); 79 PetscCall(PetscFree(islocal2)); 80 } 81 PetscCall(MatDestroy(&C)); 82 PetscCall(PetscFinalize()); 83 return 0; 84 } 85 86 /*TEST 87 88 test: 89 args: -m 7 90 91 TEST*/ 92