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