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
FormElementStiffness(PetscReal H,PetscScalar * Ke)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
main(int argc,char ** args)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