1 static char help[] = "Check to see of DM_BOUNDARY_MIRROR works in 3D for DMDA with star stencil\n";
2
3 #include "petscdmda.h"
4
5 /* Contributed by Gourav Kumbhojkar */
6
globalKMat_3d(Mat K,DMDALocalInfo info)7 PetscErrorCode globalKMat_3d(Mat K, DMDALocalInfo info)
8 {
9 MatStencil row, col[7];
10 PetscScalar vals[7];
11 PetscInt ncols;
12
13 PetscFunctionBeginUser;
14 for (PetscInt i = info.xs; i < info.xs + info.xm; i++) {
15 for (PetscInt j = info.ys; j < info.ys + info.ym; j++) {
16 for (PetscInt k = info.zs; k < info.zs + info.zm; k++) {
17 ncols = 0;
18 row.i = i;
19 row.j = j;
20 row.k = k;
21
22 col[0].i = i;
23 col[0].j = j;
24 col[0].k = k;
25 vals[ncols++] = -6.; //ncols=1
26
27 col[ncols].i = i - 1;
28 col[ncols].j = j;
29 col[ncols].k = k;
30 vals[ncols++] = 1.; //ncols=2
31
32 col[ncols].i = i + 1;
33 col[ncols].j = j;
34 col[ncols].k = k;
35 vals[ncols++] = 1.; //ncols=3
36
37 col[ncols].i = i;
38 col[ncols].j = j - 1;
39 col[ncols].k = k;
40 vals[ncols++] = 1.; //ncols=4
41
42 col[ncols].i = i;
43 col[ncols].j = j + 1;
44 col[ncols].k = k;
45 vals[ncols++] = 1.; //ncols=5
46
47 col[ncols].i = i;
48 col[ncols].j = j;
49 col[ncols].k = k + 1;
50 vals[ncols++] = 1.; //ncols=6
51
52 col[ncols].i = i;
53 col[ncols].j = j;
54 col[ncols].k = k - 1;
55 vals[ncols++] = 1.; //ncols=7
56
57 PetscCall(MatSetValuesStencil(K, 1, &row, ncols, col, vals, ADD_VALUES));
58 }
59 }
60 }
61 PetscCall(MatAssemblyBegin(K, MAT_FINAL_ASSEMBLY));
62 PetscCall(MatAssemblyEnd(K, MAT_FINAL_ASSEMBLY));
63 PetscFunctionReturn(PETSC_SUCCESS);
64 }
65
globalKMat_2d(Mat K,DMDALocalInfo info)66 PetscErrorCode globalKMat_2d(Mat K, DMDALocalInfo info)
67 {
68 MatStencil row, col[5];
69 PetscScalar vals[5];
70 PetscInt ncols;
71
72 PetscFunctionBeginUser;
73 for (PetscInt i = info.xs; i < info.xs + info.xm; i++) {
74 for (PetscInt j = info.ys; j < info.ys + info.ym; j++) {
75 ncols = 0;
76 row.i = i;
77 row.j = j;
78
79 col[0].i = i;
80 col[0].j = j;
81 vals[ncols++] = -4.; //ncols=1
82
83 col[ncols].i = i - 1;
84 col[ncols].j = j;
85 vals[ncols++] = 1.; //ncols=2
86
87 col[ncols].i = i;
88 col[ncols].j = j - 1;
89 vals[ncols++] = 1.; //ncols=3
90
91 col[ncols].i = i + 1;
92 col[ncols].j = j;
93 vals[ncols++] = 1.; //ncols=4
94
95 col[ncols].i = i;
96 col[ncols].j = j + 1;
97 vals[ncols++] = 1.; //ncols=5
98
99 PetscCall(MatSetValuesStencil(K, 1, &row, ncols, col, vals, ADD_VALUES));
100 }
101 }
102 PetscCall(MatAssemblyBegin(K, MAT_FINAL_ASSEMBLY));
103 PetscCall(MatAssemblyEnd(K, MAT_FINAL_ASSEMBLY));
104 PetscFunctionReturn(PETSC_SUCCESS);
105 }
106
main(int argc,char ** argv)107 int main(int argc, char **argv)
108 {
109 DM da3d, da2d;
110 DMDALocalInfo info3d, info2d;
111 Mat K3d, K2d;
112 PetscInt ne, num_pts;
113 ISLocalToGlobalMapping ltgm3d, ltgm2d;
114 Vec row2d, row3d;
115 PetscReal norm2d, norm3d;
116
117 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
118 ne = 8;
119 num_pts = ne + 1;
120
121 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_MIRROR, DM_BOUNDARY_MIRROR, DMDA_STENCIL_STAR, num_pts, num_pts + 1, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da2d));
122 PetscCall(DMSetUp(da2d));
123 PetscCall(DMDAGetLocalInfo(da2d, &info2d));
124 PetscCall(DMCreateMatrix(da2d, &K2d));
125 PetscCall(DMGetLocalToGlobalMapping(da2d, <gm2d));
126 PetscCall(ISLocalToGlobalMappingView(ltgm2d, PETSC_VIEWER_STDOUT_WORLD));
127 //PetscFinalize();
128 PetscCall(globalKMat_2d(K2d, info2d));
129 PetscCall(MatView(K2d, PETSC_VIEWER_STDOUT_WORLD));
130 PetscCall(MatCreateVecs(K2d, &row2d, NULL));
131
132 PetscCall(MatGetRowSum(K2d, row2d));
133 PetscCall(VecNorm(row2d, NORM_2, &norm2d));
134
135 PetscCheck(norm2d == 0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "2D atrix row sum should be zero");
136 PetscCall(VecDestroy(&row2d));
137 PetscCall(MatDestroy(&K2d));
138 PetscCall(DMDestroy(&da2d));
139
140 PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_MIRROR, DM_BOUNDARY_MIRROR, DM_BOUNDARY_MIRROR, DMDA_STENCIL_STAR, num_pts, num_pts + 1, num_pts + 2, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, NULL, &da3d));
141 PetscCall(DMSetUp(da3d));
142 PetscCall(DMCreateMatrix(da3d, &K3d));
143 PetscCall(DMDAGetLocalInfo(da3d, &info3d));
144 PetscCall(DMGetLocalToGlobalMapping(da3d, <gm3d));
145 PetscCall(ISLocalToGlobalMappingView(ltgm3d, PETSC_VIEWER_STDOUT_WORLD));
146 PetscCall(globalKMat_3d(K3d, info3d));
147 PetscCall(MatView(K3d, PETSC_VIEWER_STDOUT_WORLD));
148 PetscCall(MatCreateVecs(K3d, &row3d, NULL));
149 PetscCall(MatGetRowSum(K3d, row3d));
150 PetscCall(VecNorm(row3d, NORM_2, &norm3d));
151 PetscCheck(norm3d == 0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "3D atrix row sum should be zero");
152 PetscCall(VecDestroy(&row3d));
153
154 PetscCall(DMDestroy(&da3d));
155 PetscCall(MatDestroy(&K3d));
156 return PetscFinalize();
157 }
158
159 /*TEST
160
161 test:
162
163 test:
164 suffix: 2
165 nsize: 2
166
167 test:
168 suffix: 4
169 nsize: 4
170
171 test:
172 suffix: 8
173 nsize: 8
174
175 TEST*/
176