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