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