1 static char help[] = "Test VTK structured (.vts) and rectilinear (.vtr) viewer support with multi-dof DMDAs.\n\ 2 Supply the -namefields flag to test with field names.\n\n"; 3 4 #include <petscdm.h> 5 #include <petscdmda.h> 6 7 /* Helper function to name DMDA fields */ 8 PetscErrorCode NameFields(DM da, PetscInt dof) 9 { 10 PetscInt c; 11 12 PetscFunctionBeginUser; 13 for (c = 0; c < dof; ++c) { 14 char fieldname[256]; 15 PetscCall(PetscSNPrintf(fieldname, sizeof(fieldname), "field_%" PetscInt_FMT, c)); 16 PetscCall(DMDASetFieldName(da, c, fieldname)); 17 } 18 PetscFunctionReturn(PETSC_SUCCESS); 19 } 20 21 /* 22 Write 3D DMDA vector with coordinates in VTK format 23 */ 24 PetscErrorCode test_3d(const char filename[], PetscInt dof, PetscBool namefields) 25 { 26 MPI_Comm comm = MPI_COMM_WORLD; 27 const PetscInt M = 10, N = 15, P = 30, sw = 1; 28 const PetscScalar Lx = 1.0, Ly = 1.0, Lz = 1.0; 29 DM da; 30 Vec v; 31 PetscViewer view; 32 DMDALocalInfo info; 33 PetscScalar ****va; 34 PetscInt i, j, k, c; 35 36 PetscCall(DMDACreate3d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, P, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof, sw, NULL, NULL, NULL, &da)); 37 PetscCall(DMSetFromOptions(da)); 38 PetscCall(DMSetUp(da)); 39 if (namefields) PetscCall(NameFields(da, dof)); 40 41 PetscCall(DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz)); 42 PetscCall(DMDAGetLocalInfo(da, &info)); 43 PetscCall(DMCreateGlobalVector(da, &v)); 44 PetscCall(DMDAVecGetArrayDOF(da, v, &va)); 45 for (k = info.zs; k < info.zs + info.zm; k++) { 46 for (j = info.ys; j < info.ys + info.ym; j++) { 47 for (i = info.xs; i < info.xs + info.xm; i++) { 48 const PetscScalar x = (Lx * i) / M; 49 const PetscScalar y = (Ly * j) / N; 50 const PetscScalar z = (Lz * k) / P; 51 for (c = 0; c < dof; ++c) va[k][j][i][c] = PetscPowScalarInt(x - 0.5 * Lx, 2) + PetscPowScalarInt(y - 0.5 * Ly, 2) + PetscPowScalarInt(z - 0.5 * Lz, 2) + 10.0 * c; 52 } 53 } 54 } 55 PetscCall(DMDAVecRestoreArrayDOF(da, v, &va)); 56 PetscCall(PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view)); 57 PetscCall(VecView(v, view)); 58 PetscCall(PetscViewerDestroy(&view)); 59 PetscCall(VecDestroy(&v)); 60 PetscCall(DMDestroy(&da)); 61 return PETSC_SUCCESS; 62 } 63 64 /* 65 Write 2D DMDA vector with coordinates in VTK format 66 */ 67 PetscErrorCode test_2d(const char filename[], PetscInt dof, PetscBool namefields) 68 { 69 MPI_Comm comm = MPI_COMM_WORLD; 70 const PetscInt M = 10, N = 20, sw = 1; 71 const PetscScalar Lx = 1.0, Ly = 1.0, Lz = 1.0; 72 DM da; 73 Vec v; 74 PetscViewer view; 75 DMDALocalInfo info; 76 PetscScalar ***va; 77 PetscInt i, j, c; 78 79 PetscCall(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, PETSC_DECIDE, PETSC_DECIDE, dof, sw, NULL, NULL, &da)); 80 PetscCall(DMSetFromOptions(da)); 81 PetscCall(DMSetUp(da)); 82 if (namefields) PetscCall(NameFields(da, dof)); 83 PetscCall(DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz)); 84 PetscCall(DMDAGetLocalInfo(da, &info)); 85 PetscCall(DMCreateGlobalVector(da, &v)); 86 PetscCall(DMDAVecGetArrayDOF(da, v, &va)); 87 for (j = info.ys; j < info.ys + info.ym; j++) { 88 for (i = info.xs; i < info.xs + info.xm; i++) { 89 const PetscScalar x = (Lx * i) / M; 90 const PetscScalar y = (Ly * j) / N; 91 for (c = 0; c < dof; ++c) va[j][i][c] = PetscPowScalarInt(x - 0.5 * Lx, 2) + PetscPowScalarInt(y - 0.5 * Ly, 2) + 10.0 * c; 92 } 93 } 94 PetscCall(DMDAVecRestoreArrayDOF(da, v, &va)); 95 PetscCall(PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view)); 96 PetscCall(VecView(v, view)); 97 PetscCall(PetscViewerDestroy(&view)); 98 PetscCall(VecDestroy(&v)); 99 PetscCall(DMDestroy(&da)); 100 return PETSC_SUCCESS; 101 } 102 103 /* 104 Write a scalar and a vector field from two compatible 3d DMDAs 105 */ 106 PetscErrorCode test_3d_compat(const char filename[], PetscInt dof, PetscBool namefields) 107 { 108 MPI_Comm comm = MPI_COMM_WORLD; 109 const PetscInt M = 10, N = 15, P = 30, sw = 1; 110 const PetscScalar Lx = 1.0, Ly = 1.0, Lz = 1.0; 111 DM da, daVector; 112 Vec v, vVector; 113 PetscViewer view; 114 DMDALocalInfo info; 115 PetscScalar ***va, ****vVectora; 116 PetscInt i, j, k, c; 117 118 PetscCall(DMDACreate3d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, P, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, /* dof:*/ 1, sw, NULL, NULL, NULL, &da)); 119 PetscCall(DMSetFromOptions(da)); 120 PetscCall(DMSetUp(da)); 121 if (namefields) PetscCall(NameFields(da, 1)); 122 123 PetscCall(DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz)); 124 PetscCall(DMDAGetLocalInfo(da, &info)); 125 PetscCall(DMDACreateCompatibleDMDA(da, dof, &daVector)); 126 if (namefields) PetscCall(NameFields(daVector, dof)); 127 PetscCall(DMCreateGlobalVector(da, &v)); 128 PetscCall(DMCreateGlobalVector(daVector, &vVector)); 129 PetscCall(DMDAVecGetArray(da, v, &va)); 130 PetscCall(DMDAVecGetArrayDOF(daVector, vVector, &vVectora)); 131 for (k = info.zs; k < info.zs + info.zm; k++) { 132 for (j = info.ys; j < info.ys + info.ym; j++) { 133 for (i = info.xs; i < info.xs + info.xm; i++) { 134 const PetscScalar x = (Lx * i) / M; 135 const PetscScalar y = (Ly * j) / N; 136 const PetscScalar z = (Lz * k) / P; 137 va[k][j][i] = PetscPowScalarInt(x - 0.5 * Lx, 2) + PetscPowScalarInt(y - 0.5 * Ly, 2) + PetscPowScalarInt(z - 0.5 * Lz, 2); 138 for (c = 0; c < dof; ++c) vVectora[k][j][i][c] = PetscPowScalarInt(x - 0.5 * Lx, 2) + PetscPowScalarInt(y - 0.5 * Ly, 2) + PetscPowScalarInt(z - 0.5 * Lz, 2) + 10.0 * c; 139 } 140 } 141 } 142 PetscCall(DMDAVecRestoreArray(da, v, &va)); 143 PetscCall(DMDAVecRestoreArrayDOF(da, v, &vVectora)); 144 PetscCall(PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view)); 145 PetscCall(VecView(v, view)); 146 PetscCall(VecView(vVector, view)); 147 PetscCall(PetscViewerDestroy(&view)); 148 PetscCall(VecDestroy(&v)); 149 PetscCall(VecDestroy(&vVector)); 150 PetscCall(DMDestroy(&da)); 151 PetscCall(DMDestroy(&daVector)); 152 return PETSC_SUCCESS; 153 } 154 155 /* 156 Write a scalar and a vector field from two compatible 2d DMDAs 157 */ 158 PetscErrorCode test_2d_compat(const char filename[], PetscInt dof, PetscBool namefields) 159 { 160 MPI_Comm comm = MPI_COMM_WORLD; 161 const PetscInt M = 10, N = 20, sw = 1; 162 const PetscScalar Lx = 1.0, Ly = 1.0, Lz = 1.0; 163 DM da, daVector; 164 Vec v, vVector; 165 PetscViewer view; 166 DMDALocalInfo info; 167 PetscScalar **va, ***vVectora; 168 PetscInt i, j, c; 169 170 PetscCall(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, PETSC_DECIDE, PETSC_DECIDE, /* dof:*/ 1, sw, NULL, NULL, &da)); 171 PetscCall(DMSetFromOptions(da)); 172 PetscCall(DMSetUp(da)); 173 if (namefields) PetscCall(NameFields(da, 1)); 174 PetscCall(DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz)); 175 PetscCall(DMDACreateCompatibleDMDA(da, dof, &daVector)); 176 if (namefields) PetscCall(NameFields(daVector, dof)); 177 PetscCall(DMDAGetLocalInfo(da, &info)); 178 PetscCall(DMCreateGlobalVector(da, &v)); 179 PetscCall(DMCreateGlobalVector(daVector, &vVector)); 180 PetscCall(DMDAVecGetArray(da, v, &va)); 181 PetscCall(DMDAVecGetArrayDOF(daVector, vVector, &vVectora)); 182 for (j = info.ys; j < info.ys + info.ym; j++) { 183 for (i = info.xs; i < info.xs + info.xm; i++) { 184 const PetscScalar x = (Lx * i) / M; 185 const PetscScalar y = (Ly * j) / N; 186 va[j][i] = PetscPowScalarInt(x - 0.5 * Lx, 2) + PetscPowScalarInt(y - 0.5 * Ly, 2); 187 for (c = 0; c < dof; ++c) vVectora[j][i][c] = PetscPowScalarInt(x - 0.5 * Lx, 2) + PetscPowScalarInt(y - 0.5 * Ly, 2) + 10.0 * c; 188 } 189 } 190 PetscCall(DMDAVecRestoreArray(da, v, &va)); 191 PetscCall(DMDAVecRestoreArrayDOF(daVector, vVector, &vVectora)); 192 PetscCall(PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view)); 193 PetscCall(VecView(v, view)); 194 PetscCall(VecView(vVector, view)); 195 PetscCall(PetscViewerDestroy(&view)); 196 PetscCall(VecDestroy(&v)); 197 PetscCall(VecDestroy(&vVector)); 198 PetscCall(DMDestroy(&da)); 199 PetscCall(DMDestroy(&daVector)); 200 return PETSC_SUCCESS; 201 } 202 203 int main(int argc, char *argv[]) 204 { 205 PetscInt dof; 206 PetscBool namefields; 207 208 PetscFunctionBeginUser; 209 PetscCall(PetscInitialize(&argc, &argv, 0, help)); 210 dof = 2; 211 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dof", &dof, NULL)); 212 namefields = PETSC_FALSE; 213 PetscCall(PetscOptionsGetBool(NULL, NULL, "-namefields", &namefields, NULL)); 214 PetscCall(test_3d("3d.vtr", dof, namefields)); 215 PetscCall(test_2d("2d.vtr", dof, namefields)); 216 PetscCall(test_3d_compat("3d_compat.vtr", dof, namefields)); 217 PetscCall(test_2d_compat("2d_compat.vtr", dof, namefields)); 218 PetscCall(test_3d("3d.vts", dof, namefields)); 219 PetscCall(test_2d("2d.vts", dof, namefields)); 220 PetscCall(test_3d_compat("3d_compat.vts", dof, namefields)); 221 PetscCall(test_2d_compat("2d_compat.vts", dof, namefields)); 222 PetscCall(PetscFinalize()); 223 return 0; 224 } 225 226 /*TEST 227 228 build: 229 requires: !complex 230 231 test: 232 suffix: 1 233 nsize: 2 234 args: -dof 2 235 236 test: 237 suffix: 2 238 nsize: 2 239 args: -dof 2 240 241 test: 242 suffix: 3 243 nsize: 2 244 args: -dof 3 245 246 test: 247 suffix: 4 248 nsize: 1 249 args: -dof 2 -namefields 250 251 test: 252 suffix: 5 253 nsize: 2 254 args: -dof 4 -namefields 255 256 TEST*/ 257