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 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 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 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, (char *)0, 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