Actual source code: ex6.c
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";
6: #include <petscksp.h>
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: }
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
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;
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));
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));
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));
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"));
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]));
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: }
105: /*TEST
107: test:
108: args: -m 7
110: TEST*/