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