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