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 8*9371c9d4SSatish Balay PetscErrorCode FormElementStiffness(PetscReal H, PetscScalar *Ke) { 9*9371c9d4SSatish Balay Ke[0] = H / 6.0; 10*9371c9d4SSatish Balay Ke[1] = -.125 * H; 11*9371c9d4SSatish Balay Ke[2] = H / 12.0; 12*9371c9d4SSatish Balay Ke[3] = -.125 * H; 13*9371c9d4SSatish Balay Ke[4] = -.125 * H; 14*9371c9d4SSatish Balay Ke[5] = H / 6.0; 15*9371c9d4SSatish Balay Ke[6] = -.125 * H; 16*9371c9d4SSatish Balay Ke[7] = H / 12.0; 17*9371c9d4SSatish Balay Ke[8] = H / 12.0; 18*9371c9d4SSatish Balay Ke[9] = -.125 * H; 19*9371c9d4SSatish Balay Ke[10] = H / 6.0; 20*9371c9d4SSatish Balay Ke[11] = -.125 * H; 21*9371c9d4SSatish Balay Ke[12] = -.125 * H; 22*9371c9d4SSatish Balay Ke[13] = H / 12.0; 23*9371c9d4SSatish Balay Ke[14] = -.125 * H; 24*9371c9d4SSatish Balay Ke[15] = H / 6.0; 25c4762a1bSJed Brown return 0; 26c4762a1bSJed Brown } 27*9371c9d4SSatish Balay PetscErrorCode FormElementRhs(PetscReal x, PetscReal y, PetscReal H, PetscScalar *r) { 28*9371c9d4SSatish Balay r[0] = 0.; 29*9371c9d4SSatish Balay r[1] = 0.; 30*9371c9d4SSatish Balay r[2] = 0.; 31*9371c9d4SSatish Balay r[3] = 0.0; 32c4762a1bSJed Brown return 0; 33c4762a1bSJed Brown } 34c4762a1bSJed Brown 35*9371c9d4SSatish Balay int main(int argc, char **args) { 36c4762a1bSJed Brown Mat C; 37c4762a1bSJed Brown PetscInt i, m = 2, N, M, idx[4], Nsub1, Nsub2, ol = 1, x1, x2; 38c4762a1bSJed Brown PetscScalar Ke[16]; 39c4762a1bSJed Brown PetscReal h; 40c4762a1bSJed Brown IS *is1, *is2, *islocal1, *islocal2; 41c4762a1bSJed Brown PetscBool flg; 42c4762a1bSJed Brown 43327415f7SBarry Smith PetscFunctionBeginUser; 449566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 459566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL)); 46c4762a1bSJed Brown N = (m + 1) * (m + 1); /* dimension of matrix */ 47c4762a1bSJed Brown M = m * m; /* number of elements */ 48c4762a1bSJed Brown h = 1.0 / m; /* mesh width */ 49c4762a1bSJed Brown x1 = (m + 1) / 2; 50c4762a1bSJed Brown x2 = x1; 519566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL, NULL, "-x1", &x1, NULL)); 529566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL, NULL, "-x2", &x2, NULL)); 53c4762a1bSJed Brown /* create stiffness matrix */ 549566063dSJacob Faibussowitsch PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, N, N, 9, NULL, &C)); 55c4762a1bSJed Brown 56c4762a1bSJed Brown /* forms the element stiffness for the Laplacian */ 579566063dSJacob Faibussowitsch PetscCall(FormElementStiffness(h * h, Ke)); 58c4762a1bSJed Brown for (i = 0; i < M; i++) { 59c4762a1bSJed Brown /* node numbers for the four corners of element */ 60c4762a1bSJed Brown idx[0] = (m + 1) * (i / m) + (i % m); 61*9371c9d4SSatish Balay idx[1] = idx[0] + 1; 62*9371c9d4SSatish Balay idx[2] = idx[1] + m + 1; 63*9371c9d4SSatish Balay idx[3] = idx[2] - 1; 649566063dSJacob Faibussowitsch PetscCall(MatSetValues(C, 4, idx, 4, idx, Ke, ADD_VALUES)); 65c4762a1bSJed Brown } 669566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); 679566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); 68c4762a1bSJed Brown 69c4762a1bSJed Brown for (ol = 0; ol < m + 2; ++ol) { 709566063dSJacob Faibussowitsch PetscCall(PCASMCreateSubdomains2D(m + 1, m + 1, x1, x2, 1, 0, &Nsub1, &is1, &islocal1)); 719566063dSJacob Faibussowitsch PetscCall(MatIncreaseOverlap(C, Nsub1, is1, ol)); 729566063dSJacob Faibussowitsch PetscCall(PCASMCreateSubdomains2D(m + 1, m + 1, x1, x2, 1, ol, &Nsub2, &is2, &islocal2)); 73c4762a1bSJed Brown 749566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_SELF, "flg == 1 => both index sets are same\n")); 75*9371c9d4SSatish Balay if (Nsub1 != Nsub2) { PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: No of indes sets don't match\n")); } 76c4762a1bSJed Brown 77c4762a1bSJed Brown for (i = 0; i < Nsub1; ++i) { 789566063dSJacob Faibussowitsch PetscCall(ISEqual(is1[i], is2[i], &flg)); 7963a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_SELF, "i = %" PetscInt_FMT ",flg = %d \n", i, (int)flg)); 80c4762a1bSJed Brown } 819566063dSJacob Faibussowitsch for (i = 0; i < Nsub1; ++i) PetscCall(ISDestroy(&is1[i])); 829566063dSJacob Faibussowitsch for (i = 0; i < Nsub2; ++i) PetscCall(ISDestroy(&is2[i])); 839566063dSJacob Faibussowitsch for (i = 0; i < Nsub1; ++i) PetscCall(ISDestroy(&islocal1[i])); 849566063dSJacob Faibussowitsch for (i = 0; i < Nsub2; ++i) PetscCall(ISDestroy(&islocal2[i])); 85c4762a1bSJed Brown 869566063dSJacob Faibussowitsch PetscCall(PetscFree(is1)); 879566063dSJacob Faibussowitsch PetscCall(PetscFree(is2)); 889566063dSJacob Faibussowitsch PetscCall(PetscFree(islocal1)); 899566063dSJacob Faibussowitsch PetscCall(PetscFree(islocal2)); 90c4762a1bSJed Brown } 919566063dSJacob Faibussowitsch PetscCall(MatDestroy(&C)); 929566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 93b122ec5aSJacob Faibussowitsch return 0; 94c4762a1bSJed Brown } 95c4762a1bSJed Brown 96c4762a1bSJed Brown /*TEST 97c4762a1bSJed Brown 98c4762a1bSJed Brown test: 99c4762a1bSJed Brown args: -m 7 100c4762a1bSJed Brown 101c4762a1bSJed Brown TEST*/ 102