Actual source code: noflux_check.c
1: static char help[] = "Check to see of DM_BOUNDARY_MIRROR works in 3D for DMDA with star stencil\n";
3: #include "petscdmda.h"
5: /* Contributed by Gourav Kumbhojkar */
7: PetscErrorCode globalKMat_3d(Mat K, DMDALocalInfo info)
8: {
9: MatStencil row, col[7];
10: PetscScalar vals[7];
11: PetscInt ncols;
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;
22: col[0].i = i;
23: col[0].j = j;
24: col[0].k = k;
25: vals[ncols++] = -6.; //ncols=1
27: col[ncols].i = i - 1;
28: col[ncols].j = j;
29: col[ncols].k = k;
30: vals[ncols++] = 1.; //ncols=2
32: col[ncols].i = i + 1;
33: col[ncols].j = j;
34: col[ncols].k = k;
35: vals[ncols++] = 1.; //ncols=3
37: col[ncols].i = i;
38: col[ncols].j = j - 1;
39: col[ncols].k = k;
40: vals[ncols++] = 1.; //ncols=4
42: col[ncols].i = i;
43: col[ncols].j = j + 1;
44: col[ncols].k = k;
45: vals[ncols++] = 1.; //ncols=5
47: col[ncols].i = i;
48: col[ncols].j = j;
49: col[ncols].k = k + 1;
50: vals[ncols++] = 1.; //ncols=6
52: col[ncols].i = i;
53: col[ncols].j = j;
54: col[ncols].k = k - 1;
55: vals[ncols++] = 1.; //ncols=7
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: }
66: PetscErrorCode globalKMat_2d(Mat K, DMDALocalInfo info)
67: {
68: MatStencil row, col[5];
69: PetscScalar vals[5];
70: PetscInt ncols;
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;
79: col[0].i = i;
80: col[0].j = j;
81: vals[ncols++] = -4.; //ncols=1
83: col[ncols].i = i - 1;
84: col[ncols].j = j;
85: vals[ncols++] = 1.; //ncols=2
87: col[ncols].i = i;
88: col[ncols].j = j - 1;
89: vals[ncols++] = 1.; //ncols=3
91: col[ncols].i = i + 1;
92: col[ncols].j = j;
93: vals[ncols++] = 1.; //ncols=4
95: col[ncols].i = i;
96: col[ncols].j = j + 1;
97: vals[ncols++] = 1.; //ncols=5
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: }
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;
117: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
118: ne = 8;
119: num_pts = ne + 1;
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));
132: PetscCall(MatGetRowSum(K2d, row2d));
133: PetscCall(VecNorm(row2d, NORM_2, &norm2d));
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));
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));
154: PetscCall(DMDestroy(&da3d));
155: PetscCall(MatDestroy(&K3d));
156: return PetscFinalize();
157: }
159: /*TEST
161: test:
163: test:
164: suffix: 2
165: nsize: 2
167: test:
168: suffix: 4
169: nsize: 4
171: test:
172: suffix: 8
173: nsize: 8
175: TEST*/