1c4762a1bSJed Brown 2c4762a1bSJed Brown static char help[] = "Creates a matrix using 9 pt stencil, and uses it to test MatIncreaseOverlap (needed for additive Schwarz preconditioner). \n\ 3c4762a1bSJed Brown -m <size> : problem size\n\ 4c4762a1bSJed Brown -x1, -x2 <size> : no of subdomains in x and y directions\n\n"; 5c4762a1bSJed Brown 6c4762a1bSJed Brown #include <petscksp.h> 7c4762a1bSJed Brown 8c4762a1bSJed Brown PetscErrorCode FormElementStiffness(PetscReal H,PetscScalar *Ke) 9c4762a1bSJed Brown { 10c4762a1bSJed Brown Ke[0] = H/6.0; Ke[1] = -.125*H; Ke[2] = H/12.0; Ke[3] = -.125*H; 11c4762a1bSJed Brown Ke[4] = -.125*H; Ke[5] = H/6.0; Ke[6] = -.125*H; Ke[7] = H/12.0; 12c4762a1bSJed Brown Ke[8] = H/12.0; Ke[9] = -.125*H; Ke[10] = H/6.0; Ke[11] = -.125*H; 13c4762a1bSJed Brown Ke[12] = -.125*H; Ke[13] = H/12.0; Ke[14] = -.125*H; Ke[15] = H/6.0; 14c4762a1bSJed Brown return 0; 15c4762a1bSJed Brown } 16c4762a1bSJed Brown PetscErrorCode FormElementRhs(PetscReal x,PetscReal y,PetscReal H,PetscScalar *r) 17c4762a1bSJed Brown { 18c4762a1bSJed Brown r[0] = 0.; r[1] = 0.; r[2] = 0.; r[3] = 0.0; 19c4762a1bSJed Brown return 0; 20c4762a1bSJed Brown } 21c4762a1bSJed Brown 22c4762a1bSJed Brown int main(int argc,char **args) 23c4762a1bSJed Brown { 24c4762a1bSJed Brown Mat C; 25c4762a1bSJed Brown PetscInt i,m = 2,N,M,idx[4],Nsub1,Nsub2,ol=1,x1,x2; 26c4762a1bSJed Brown PetscScalar Ke[16]; 27c4762a1bSJed Brown PetscReal h; 28c4762a1bSJed Brown IS *is1,*is2,*islocal1,*islocal2; 29c4762a1bSJed Brown PetscBool flg; 30c4762a1bSJed Brown 31*b122ec5aSJacob Faibussowitsch CHKERRQ(PetscInitialize(&argc,&args,(char*)0,help)); 325f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 33c4762a1bSJed Brown N = (m+1)*(m+1); /* dimension of matrix */ 34c4762a1bSJed Brown M = m*m; /* number of elements */ 35c4762a1bSJed Brown h = 1.0/m; /* mesh width */ 36c4762a1bSJed Brown x1 = (m+1)/2; 37c4762a1bSJed Brown x2 = x1; 385f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-x1",&x1,NULL)); 395f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-x2",&x2,NULL)); 40c4762a1bSJed Brown /* create stiffness matrix */ 415f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,9,NULL,&C)); 42c4762a1bSJed Brown 43c4762a1bSJed Brown /* forms the element stiffness for the Laplacian */ 445f80ce2aSJacob Faibussowitsch CHKERRQ(FormElementStiffness(h*h,Ke)); 45c4762a1bSJed Brown for (i=0; i<M; i++) { 46c4762a1bSJed Brown /* node numbers for the four corners of element */ 47c4762a1bSJed Brown idx[0] = (m+1)*(i/m) + (i % m); 48c4762a1bSJed Brown idx[1] = idx[0]+1; idx[2] = idx[1] + m + 1; idx[3] = idx[2] - 1; 495f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetValues(C,4,idx,4,idx,Ke,ADD_VALUES)); 50c4762a1bSJed Brown } 515f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 525f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 53c4762a1bSJed Brown 54c4762a1bSJed Brown for (ol=0; ol<m+2; ++ol) { 55c4762a1bSJed Brown 565f80ce2aSJacob Faibussowitsch CHKERRQ(PCASMCreateSubdomains2D(m+1,m+1,x1,x2,1,0,&Nsub1,&is1,&islocal1)); 575f80ce2aSJacob Faibussowitsch CHKERRQ(MatIncreaseOverlap(C,Nsub1,is1,ol)); 585f80ce2aSJacob Faibussowitsch CHKERRQ(PCASMCreateSubdomains2D(m+1,m+1,x1,x2,1,ol,&Nsub2,&is2,&islocal2)); 59c4762a1bSJed Brown 605f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"flg == 1 => both index sets are same\n")); 61c4762a1bSJed Brown if (Nsub1 != Nsub2) { 625f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error: No of indes sets don't match\n")); 63c4762a1bSJed Brown } 64c4762a1bSJed Brown 65c4762a1bSJed Brown for (i=0; i<Nsub1; ++i) { 665f80ce2aSJacob Faibussowitsch CHKERRQ(ISEqual(is1[i],is2[i],&flg)); 675f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"i = %D,flg = %d \n",i,(int)flg)); 68c4762a1bSJed Brown 69c4762a1bSJed Brown } 705f80ce2aSJacob Faibussowitsch for (i=0; i<Nsub1; ++i) CHKERRQ(ISDestroy(&is1[i])); 715f80ce2aSJacob Faibussowitsch for (i=0; i<Nsub2; ++i) CHKERRQ(ISDestroy(&is2[i])); 725f80ce2aSJacob Faibussowitsch for (i=0; i<Nsub1; ++i) CHKERRQ(ISDestroy(&islocal1[i])); 735f80ce2aSJacob Faibussowitsch for (i=0; i<Nsub2; ++i) CHKERRQ(ISDestroy(&islocal2[i])); 74c4762a1bSJed Brown 755f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFree(is1)); 765f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFree(is2)); 775f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFree(islocal1)); 785f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFree(islocal2)); 79c4762a1bSJed Brown } 805f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&C)); 81*b122ec5aSJacob Faibussowitsch CHKERRQ(PetscFinalize()); 82*b122ec5aSJacob Faibussowitsch return 0; 83c4762a1bSJed Brown } 84c4762a1bSJed Brown 85c4762a1bSJed Brown /*TEST 86c4762a1bSJed Brown 87c4762a1bSJed Brown test: 88c4762a1bSJed Brown args: -m 7 89c4762a1bSJed Brown 90c4762a1bSJed Brown TEST*/ 91