1*c4762a1bSJed Brown 2*c4762a1bSJed Brown static char help[] = "Creates a matrix using 9 pt stencil, and uses it to test MatIncreaseOverlap (needed for additive Schwarz preconditioner). \n\ 3*c4762a1bSJed Brown -m <size> : problem size\n\ 4*c4762a1bSJed Brown -x1, -x2 <size> : no of subdomains in x and y directions\n\n"; 5*c4762a1bSJed Brown 6*c4762a1bSJed Brown #include <petscksp.h> 7*c4762a1bSJed Brown 8*c4762a1bSJed Brown PetscErrorCode FormElementStiffness(PetscReal H,PetscScalar *Ke) 9*c4762a1bSJed Brown { 10*c4762a1bSJed Brown Ke[0] = H/6.0; Ke[1] = -.125*H; Ke[2] = H/12.0; Ke[3] = -.125*H; 11*c4762a1bSJed Brown Ke[4] = -.125*H; Ke[5] = H/6.0; Ke[6] = -.125*H; Ke[7] = H/12.0; 12*c4762a1bSJed Brown Ke[8] = H/12.0; Ke[9] = -.125*H; Ke[10] = H/6.0; Ke[11] = -.125*H; 13*c4762a1bSJed Brown Ke[12] = -.125*H; Ke[13] = H/12.0; Ke[14] = -.125*H; Ke[15] = H/6.0; 14*c4762a1bSJed Brown return 0; 15*c4762a1bSJed Brown } 16*c4762a1bSJed Brown PetscErrorCode FormElementRhs(PetscReal x,PetscReal y,PetscReal H,PetscScalar *r) 17*c4762a1bSJed Brown { 18*c4762a1bSJed Brown r[0] = 0.; r[1] = 0.; r[2] = 0.; r[3] = 0.0; 19*c4762a1bSJed Brown return 0; 20*c4762a1bSJed Brown } 21*c4762a1bSJed Brown 22*c4762a1bSJed Brown int main(int argc,char **args) 23*c4762a1bSJed Brown { 24*c4762a1bSJed Brown Mat C; 25*c4762a1bSJed Brown PetscErrorCode ierr; 26*c4762a1bSJed Brown PetscInt i,m = 2,N,M,idx[4],Nsub1,Nsub2,ol=1,x1,x2; 27*c4762a1bSJed Brown PetscScalar Ke[16]; 28*c4762a1bSJed Brown PetscReal h; 29*c4762a1bSJed Brown IS *is1,*is2,*islocal1,*islocal2; 30*c4762a1bSJed Brown PetscBool flg; 31*c4762a1bSJed Brown 32*c4762a1bSJed Brown ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 33*c4762a1bSJed Brown ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr); 34*c4762a1bSJed Brown N = (m+1)*(m+1); /* dimension of matrix */ 35*c4762a1bSJed Brown M = m*m; /* number of elements */ 36*c4762a1bSJed Brown h = 1.0/m; /* mesh width */ 37*c4762a1bSJed Brown x1 = (m+1)/2; 38*c4762a1bSJed Brown x2 = x1; 39*c4762a1bSJed Brown ierr = PetscOptionsGetInt(NULL,NULL,"-x1",&x1,NULL);CHKERRQ(ierr); 40*c4762a1bSJed Brown ierr = PetscOptionsGetInt(NULL,NULL,"-x2",&x2,NULL);CHKERRQ(ierr); 41*c4762a1bSJed Brown /* create stiffness matrix */ 42*c4762a1bSJed Brown ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,9,NULL,&C);CHKERRQ(ierr); 43*c4762a1bSJed Brown 44*c4762a1bSJed Brown /* forms the element stiffness for the Laplacian */ 45*c4762a1bSJed Brown ierr = FormElementStiffness(h*h,Ke);CHKERRQ(ierr); 46*c4762a1bSJed Brown for (i=0; i<M; i++) { 47*c4762a1bSJed Brown /* node numbers for the four corners of element */ 48*c4762a1bSJed Brown idx[0] = (m+1)*(i/m) + (i % m); 49*c4762a1bSJed Brown idx[1] = idx[0]+1; idx[2] = idx[1] + m + 1; idx[3] = idx[2] - 1; 50*c4762a1bSJed Brown ierr = MatSetValues(C,4,idx,4,idx,Ke,ADD_VALUES);CHKERRQ(ierr); 51*c4762a1bSJed Brown } 52*c4762a1bSJed Brown ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 53*c4762a1bSJed Brown ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 54*c4762a1bSJed Brown 55*c4762a1bSJed Brown for (ol=0; ol<m+2; ++ol) { 56*c4762a1bSJed Brown 57*c4762a1bSJed Brown ierr = PCASMCreateSubdomains2D(m+1,m+1,x1,x2,1,0,&Nsub1,&is1,&islocal1);CHKERRQ(ierr); 58*c4762a1bSJed Brown ierr = MatIncreaseOverlap(C,Nsub1,is1,ol);CHKERRQ(ierr); 59*c4762a1bSJed Brown ierr = PCASMCreateSubdomains2D(m+1,m+1,x1,x2,1,ol,&Nsub2,&is2,&islocal2);CHKERRQ(ierr); 60*c4762a1bSJed Brown 61*c4762a1bSJed Brown ierr = PetscPrintf(PETSC_COMM_SELF,"flg == 1 => both index sets are same\n");CHKERRQ(ierr); 62*c4762a1bSJed Brown if (Nsub1 != Nsub2) { 63*c4762a1bSJed Brown ierr = PetscPrintf(PETSC_COMM_SELF,"Error: No of indes sets don't match\n");CHKERRQ(ierr); 64*c4762a1bSJed Brown } 65*c4762a1bSJed Brown 66*c4762a1bSJed Brown for (i=0; i<Nsub1; ++i) { 67*c4762a1bSJed Brown ierr = ISEqual(is1[i],is2[i],&flg);CHKERRQ(ierr); 68*c4762a1bSJed Brown ierr = PetscPrintf(PETSC_COMM_SELF,"i = %D,flg = %d \n",i,(int)flg);CHKERRQ(ierr); 69*c4762a1bSJed Brown 70*c4762a1bSJed Brown } 71*c4762a1bSJed Brown for (i=0; i<Nsub1; ++i) {ierr = ISDestroy(&is1[i]);CHKERRQ(ierr);} 72*c4762a1bSJed Brown for (i=0; i<Nsub2; ++i) {ierr = ISDestroy(&is2[i]);CHKERRQ(ierr);} 73*c4762a1bSJed Brown for (i=0; i<Nsub1; ++i) {ierr = ISDestroy(&islocal1[i]);CHKERRQ(ierr);} 74*c4762a1bSJed Brown for (i=0; i<Nsub2; ++i) {ierr = ISDestroy(&islocal2[i]);CHKERRQ(ierr);} 75*c4762a1bSJed Brown 76*c4762a1bSJed Brown 77*c4762a1bSJed Brown ierr = PetscFree(is1);CHKERRQ(ierr); 78*c4762a1bSJed Brown ierr = PetscFree(is2);CHKERRQ(ierr); 79*c4762a1bSJed Brown ierr = PetscFree(islocal1);CHKERRQ(ierr); 80*c4762a1bSJed Brown ierr = PetscFree(islocal2);CHKERRQ(ierr); 81*c4762a1bSJed Brown } 82*c4762a1bSJed Brown ierr = MatDestroy(&C);CHKERRQ(ierr); 83*c4762a1bSJed Brown ierr = PetscFinalize(); 84*c4762a1bSJed Brown return ierr; 85*c4762a1bSJed Brown } 86*c4762a1bSJed Brown 87*c4762a1bSJed Brown 88*c4762a1bSJed Brown 89*c4762a1bSJed Brown /*TEST 90*c4762a1bSJed Brown 91*c4762a1bSJed Brown test: 92*c4762a1bSJed Brown args: -m 7 93*c4762a1bSJed Brown 94*c4762a1bSJed Brown TEST*/ 95