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: 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: }
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;
47: PetscInitialize(&argc, &args, (char *)0, help);
48: 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: PetscOptionsGetInt(NULL, NULL, "-x1", &x1, NULL);
55: PetscOptionsGetInt(NULL, NULL, "-x2", &x2, NULL);
56: /* create stiffness matrix */
57: MatCreateSeqAIJ(PETSC_COMM_SELF, N, N, 9, NULL, &C);
59: /* forms the element stiffness for the Laplacian */
60: 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: MatSetValues(C, 4, idx, 4, idx, Ke, ADD_VALUES);
68: }
69: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
70: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
72: for (ol = 0; ol < m + 2; ++ol) {
73: PCASMCreateSubdomains2D(m + 1, m + 1, x1, x2, 1, 0, &Nsub1, &is1, &islocal1);
74: MatIncreaseOverlap(C, Nsub1, is1, ol);
75: PCASMCreateSubdomains2D(m + 1, m + 1, x1, x2, 1, ol, &Nsub2, &is2, &islocal2);
77: PetscPrintf(PETSC_COMM_SELF, "flg == 1 => both index sets are same\n");
78: if (Nsub1 != Nsub2) PetscPrintf(PETSC_COMM_SELF, "Error: No of index sets don't match\n");
80: for (i = 0; i < Nsub1; ++i) {
81: ISEqual(is1[i], is2[i], &flg);
82: PetscPrintf(PETSC_COMM_SELF, "i = %" PetscInt_FMT ",flg = %d \n", i, (int)flg);
83: }
84: for (i = 0; i < Nsub1; ++i) ISDestroy(&is1[i]);
85: for (i = 0; i < Nsub2; ++i) ISDestroy(&is2[i]);
86: for (i = 0; i < Nsub1; ++i) ISDestroy(&islocal1[i]);
87: for (i = 0; i < Nsub2; ++i) ISDestroy(&islocal2[i]);
89: PetscFree(is1);
90: PetscFree(is2);
91: PetscFree(islocal1);
92: PetscFree(islocal2);
93: }
94: MatDestroy(&C);
95: PetscFinalize();
96: return 0;
97: }
99: /*TEST
101: test:
102: args: -m 7
104: TEST*/