1 #include <petsc/private/dmdaimpl.h> 2 /* 3 Note that the API for using PETSCVIEWERVTK is totally wrong since its use requires 4 including the private vtkvimpl.h file. The code should be refactored. 5 */ 6 #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h> 7 8 /* Helper function which determines if any DMDA fields are named. This is used 9 as a proxy for the user's intention to use DMDA fields as distinct 10 scalar-valued fields as opposed to a single vector-valued field */ 11 static PetscErrorCode DMDAGetFieldsNamed(DM da, PetscBool *fieldsnamed) 12 { 13 PetscInt f, bs; 14 15 PetscFunctionBegin; 16 *fieldsnamed = PETSC_FALSE; 17 PetscCall(DMDAGetDof(da, &bs)); 18 for (f = 0; f < bs; ++f) { 19 const char *fieldname; 20 PetscCall(DMDAGetFieldName(da, f, &fieldname)); 21 if (fieldname) { 22 *fieldsnamed = PETSC_TRUE; 23 break; 24 } 25 } 26 PetscFunctionReturn(PETSC_SUCCESS); 27 } 28 29 static PetscErrorCode DMDAVTKWriteAll_VTS(DM da, PetscViewer viewer) 30 { 31 const char *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian"; 32 #if defined(PETSC_USE_REAL_SINGLE) 33 const char precision[] = "Float32"; 34 #elif defined(PETSC_USE_REAL_DOUBLE) 35 const char precision[] = "Float64"; 36 #else 37 const char precision[] = "UnknownPrecision"; 38 #endif 39 MPI_Comm comm; 40 Vec Coords; 41 PetscViewer_VTK *vtk = (PetscViewer_VTK *)viewer->data; 42 PetscViewerVTKObjectLink link; 43 FILE *fp; 44 PetscMPIInt rank, size, tag; 45 DMDALocalInfo info; 46 PetscInt dim, mx, my, mz, cdim, bs, maxnnodes, maxbs, i, j, k, r; 47 PetscInt64 boffset; 48 PetscInt rloc[6], (*grloc)[6] = NULL; 49 PetscScalar *array, *array2; 50 51 PetscFunctionBegin; 52 PetscCall(PetscObjectGetComm((PetscObject)da, &comm)); 53 PetscCheck(!PetscDefined(USE_COMPLEX), PETSC_COMM_SELF, PETSC_ERR_SUP, "Complex values not supported"); 54 PetscCallMPI(MPI_Comm_size(comm, &size)); 55 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 56 PetscCall(DMDAGetInfo(da, &dim, &mx, &my, &mz, NULL, NULL, NULL, &bs, NULL, NULL, NULL, NULL, NULL)); 57 PetscCall(DMDAGetLocalInfo(da, &info)); 58 PetscCall(DMGetCoordinates(da, &Coords)); 59 if (Coords) { 60 PetscInt csize; 61 PetscCall(VecGetSize(Coords, &csize)); 62 PetscCheck(csize % (mx * my * mz) == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Coordinate vector size mismatch"); 63 cdim = csize / (mx * my * mz); 64 PetscCheck(cdim >= dim && cdim <= 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Coordinate vector size mismatch"); 65 } else { 66 cdim = dim; 67 } 68 69 PetscCall(PetscFOpen(comm, vtk->filename, "wb", &fp)); 70 PetscCall(PetscFPrintf(comm, fp, "<?xml version=\"1.0\"?>\n")); 71 PetscCall(PetscFPrintf(comm, fp, "<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"%s\" header_type=\"UInt64\">\n", byte_order)); 72 PetscCall(PetscFPrintf(comm, fp, " <StructuredGrid WholeExtent=\"%d %" PetscInt_FMT " %d %" PetscInt_FMT " %d %" PetscInt_FMT "\">\n", 0, mx - 1, 0, my - 1, 0, mz - 1)); 73 74 if (rank == 0) PetscCall(PetscMalloc1(size * 6, &grloc)); 75 rloc[0] = info.xs; 76 rloc[1] = info.xm; 77 rloc[2] = info.ys; 78 rloc[3] = info.ym; 79 rloc[4] = info.zs; 80 rloc[5] = info.zm; 81 PetscCallMPI(MPI_Gather(rloc, 6, MPIU_INT, &grloc[0][0], 6, MPIU_INT, 0, comm)); 82 83 /* Write XML header */ 84 maxnnodes = 0; /* Used for the temporary array size on rank 0 */ 85 maxbs = 0; /* Used for the temporary array size on rank 0 */ 86 boffset = 0; /* Offset into binary file */ 87 for (r = 0; r < size; r++) { 88 PetscInt xs = -1, xm = -1, ys = -1, ym = -1, zs = -1, zm = -1, nnodes = 0; 89 if (rank == 0) { 90 xs = grloc[r][0]; 91 xm = grloc[r][1]; 92 ys = grloc[r][2]; 93 ym = grloc[r][3]; 94 zs = grloc[r][4]; 95 zm = grloc[r][5]; 96 nnodes = xm * ym * zm; 97 } 98 maxnnodes = PetscMax(maxnnodes, nnodes); 99 PetscCall(PetscFPrintf(comm, fp, " <Piece Extent=\"%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\">\n", xs, xs + xm - 1, ys, ys + ym - 1, zs, zs + zm - 1)); 100 PetscCall(PetscFPrintf(comm, fp, " <Points>\n")); 101 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, boffset)); 102 boffset += 3 * nnodes * sizeof(PetscScalar) + sizeof(PetscInt64); 103 PetscCall(PetscFPrintf(comm, fp, " </Points>\n")); 104 105 PetscCall(PetscFPrintf(comm, fp, " <PointData Scalars=\"ScalarPointData\">\n")); 106 for (link = vtk->link; link; link = link->next) { 107 Vec X = (Vec)link->vec; 108 PetscInt bs, f; 109 DM daCurr; 110 PetscBool fieldsnamed; 111 const char *vecname = "Unnamed"; 112 113 PetscCall(VecGetDM(X, &daCurr)); 114 PetscCall(DMDAGetInfo(daCurr, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bs, NULL, NULL, NULL, NULL, NULL)); 115 maxbs = PetscMax(maxbs, bs); 116 117 if (((PetscObject)X)->name || link != vtk->link) PetscCall(PetscObjectGetName((PetscObject)X, &vecname)); 118 119 /* If any fields are named, add scalar fields. Otherwise, add a vector field */ 120 PetscCall(DMDAGetFieldsNamed(daCurr, &fieldsnamed)); 121 if (fieldsnamed) { 122 for (f = 0; f < bs; f++) { 123 char buf[256]; 124 const char *fieldname; 125 PetscCall(DMDAGetFieldName(daCurr, f, &fieldname)); 126 if (!fieldname) { 127 PetscCall(PetscSNPrintf(buf, sizeof(buf), "%" PetscInt_FMT, f)); 128 fieldname = buf; 129 } 130 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s.%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset)); 131 boffset += nnodes * sizeof(PetscScalar) + sizeof(PetscInt64); 132 } 133 } else { 134 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"%" PetscInt_FMT "\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, bs, boffset)); 135 boffset += bs * nnodes * sizeof(PetscScalar) + sizeof(PetscInt64); 136 } 137 } 138 PetscCall(PetscFPrintf(comm, fp, " </PointData>\n")); 139 PetscCall(PetscFPrintf(comm, fp, " </Piece>\n")); 140 } 141 PetscCall(PetscFPrintf(comm, fp, " </StructuredGrid>\n")); 142 PetscCall(PetscFPrintf(comm, fp, " <AppendedData encoding=\"raw\">\n")); 143 PetscCall(PetscFPrintf(comm, fp, "_")); 144 145 /* Now write the arrays. */ 146 tag = ((PetscObject)viewer)->tag; 147 PetscCall(PetscMalloc2(maxnnodes * PetscMax(3, maxbs), &array, maxnnodes * PetscMax(3, maxbs), &array2)); 148 for (r = 0; r < size; r++) { 149 MPI_Status status; 150 PetscInt xs = -1, xm = -1, ys = -1, ym = -1, zs = -1, zm = -1, nnodes = 0; 151 if (rank == 0) { 152 xs = grloc[r][0]; 153 xm = grloc[r][1]; 154 ys = grloc[r][2]; 155 ym = grloc[r][3]; 156 zs = grloc[r][4]; 157 zm = grloc[r][5]; 158 nnodes = xm * ym * zm; 159 } else if (r == rank) { 160 nnodes = info.xm * info.ym * info.zm; 161 } 162 163 /* Write the coordinates */ 164 if (Coords) { 165 const PetscScalar *coords; 166 PetscCall(VecGetArrayRead(Coords, &coords)); 167 if (rank == 0) { 168 if (r) { 169 PetscMPIInt nn; 170 PetscCallMPI(MPI_Recv(array, nnodes * cdim, MPIU_SCALAR, r, tag, comm, &status)); 171 PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &nn)); 172 PetscCheck(nn == nnodes * cdim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Array size mismatch"); 173 } else PetscCall(PetscArraycpy(array, coords, nnodes * cdim)); 174 /* Transpose coordinates to VTK (C-style) ordering */ 175 for (k = 0; k < zm; k++) { 176 for (j = 0; j < ym; j++) { 177 for (i = 0; i < xm; i++) { 178 PetscInt Iloc = i + xm * (j + ym * k); 179 array2[Iloc * 3 + 0] = array[Iloc * cdim + 0]; 180 array2[Iloc * 3 + 1] = cdim > 1 ? array[Iloc * cdim + 1] : 0.0; 181 array2[Iloc * 3 + 2] = cdim > 2 ? array[Iloc * cdim + 2] : 0.0; 182 } 183 } 184 } 185 } else if (r == rank) { 186 PetscCallMPI(MPI_Send((void *)coords, nnodes * cdim, MPIU_SCALAR, 0, tag, comm)); 187 } 188 PetscCall(VecRestoreArrayRead(Coords, &coords)); 189 } else { /* Fabricate some coordinates using grid index */ 190 for (k = 0; k < zm; k++) { 191 for (j = 0; j < ym; j++) { 192 for (i = 0; i < xm; i++) { 193 PetscInt Iloc = i + xm * (j + ym * k); 194 array2[Iloc * 3 + 0] = xs + i; 195 array2[Iloc * 3 + 1] = ys + j; 196 array2[Iloc * 3 + 2] = zs + k; 197 } 198 } 199 } 200 } 201 PetscCall(PetscViewerVTKFWrite(viewer, fp, array2, nnodes * 3, MPIU_SCALAR)); 202 203 /* Write each of the objects queued up for this file */ 204 for (link = vtk->link; link; link = link->next) { 205 Vec X = (Vec)link->vec; 206 const PetscScalar *x; 207 PetscInt bs, f; 208 DM daCurr; 209 PetscBool fieldsnamed; 210 PetscCall(VecGetDM(X, &daCurr)); 211 PetscCall(DMDAGetInfo(daCurr, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bs, NULL, NULL, NULL, NULL, NULL)); 212 PetscCall(VecGetArrayRead(X, &x)); 213 if (rank == 0) { 214 if (r) { 215 PetscMPIInt nn; 216 PetscCallMPI(MPI_Recv(array, nnodes * bs, MPIU_SCALAR, r, tag, comm, &status)); 217 PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &nn)); 218 PetscCheck(nn == nnodes * bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Array size mismatch receiving from rank %" PetscInt_FMT, r); 219 } else PetscCall(PetscArraycpy(array, x, nnodes * bs)); 220 221 /* If any fields are named, add scalar fields. Otherwise, add a vector field */ 222 PetscCall(DMDAGetFieldsNamed(daCurr, &fieldsnamed)); 223 if (fieldsnamed) { 224 for (f = 0; f < bs; f++) { 225 /* Extract and transpose the f'th field */ 226 for (k = 0; k < zm; k++) { 227 for (j = 0; j < ym; j++) { 228 for (i = 0; i < xm; i++) { 229 PetscInt Iloc = i + xm * (j + ym * k); 230 array2[Iloc] = array[Iloc * bs + f]; 231 } 232 } 233 } 234 PetscCall(PetscViewerVTKFWrite(viewer, fp, array2, nnodes, MPIU_SCALAR)); 235 } 236 } else PetscCall(PetscViewerVTKFWrite(viewer, fp, array, bs * nnodes, MPIU_SCALAR)); 237 } else if (r == rank) PetscCallMPI(MPI_Send((void *)x, nnodes * bs, MPIU_SCALAR, 0, tag, comm)); 238 PetscCall(VecRestoreArrayRead(X, &x)); 239 } 240 } 241 PetscCall(PetscFree2(array, array2)); 242 PetscCall(PetscFree(grloc)); 243 244 PetscCall(PetscFPrintf(comm, fp, "\n </AppendedData>\n")); 245 PetscCall(PetscFPrintf(comm, fp, "</VTKFile>\n")); 246 PetscCall(PetscFClose(comm, fp)); 247 PetscFunctionReturn(PETSC_SUCCESS); 248 } 249 250 static PetscErrorCode DMDAVTKWriteAll_VTR(DM da, PetscViewer viewer) 251 { 252 const char *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian"; 253 #if defined(PETSC_USE_REAL_SINGLE) 254 const char precision[] = "Float32"; 255 #elif defined(PETSC_USE_REAL_DOUBLE) 256 const char precision[] = "Float64"; 257 #else 258 const char precision[] = "UnknownPrecision"; 259 #endif 260 MPI_Comm comm; 261 PetscViewer_VTK *vtk = (PetscViewer_VTK *)viewer->data; 262 PetscViewerVTKObjectLink link; 263 FILE *fp; 264 PetscMPIInt rank, size, tag; 265 DMDALocalInfo info; 266 PetscInt dim, mx, my, mz, maxnnodes, maxbs, i, j, k, r; 267 PetscInt64 boffset; 268 PetscInt rloc[6], (*grloc)[6] = NULL; 269 PetscScalar *array, *array2; 270 271 PetscFunctionBegin; 272 PetscCall(PetscObjectGetComm((PetscObject)da, &comm)); 273 PetscCheck(!PetscDefined(USE_COMPLEX), PETSC_COMM_SELF, PETSC_ERR_SUP, "Complex values not supported"); 274 PetscCallMPI(MPI_Comm_size(comm, &size)); 275 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 276 PetscCall(DMDAGetInfo(da, &dim, &mx, &my, &mz, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 277 PetscCall(DMDAGetLocalInfo(da, &info)); 278 PetscCall(PetscFOpen(comm, vtk->filename, "wb", &fp)); 279 PetscCall(PetscFPrintf(comm, fp, "<?xml version=\"1.0\"?>\n")); 280 PetscCall(PetscFPrintf(comm, fp, "<VTKFile type=\"RectilinearGrid\" version=\"0.1\" byte_order=\"%s\" header_type=\"UInt64\">\n", byte_order)); 281 PetscCall(PetscFPrintf(comm, fp, " <RectilinearGrid WholeExtent=\"%d %" PetscInt_FMT " %d %" PetscInt_FMT " %d %" PetscInt_FMT "\">\n", 0, mx - 1, 0, my - 1, 0, mz - 1)); 282 283 if (rank == 0) PetscCall(PetscMalloc1(size * 6, &grloc)); 284 rloc[0] = info.xs; 285 rloc[1] = info.xm; 286 rloc[2] = info.ys; 287 rloc[3] = info.ym; 288 rloc[4] = info.zs; 289 rloc[5] = info.zm; 290 PetscCallMPI(MPI_Gather(rloc, 6, MPIU_INT, &grloc[0][0], 6, MPIU_INT, 0, comm)); 291 292 /* Write XML header */ 293 maxnnodes = 0; /* Used for the temporary array size on rank 0 */ 294 maxbs = 0; /* Used for the temporary array size on rank 0 */ 295 boffset = 0; /* Offset into binary file */ 296 for (r = 0; r < size; r++) { 297 PetscInt xs = -1, xm = -1, ys = -1, ym = -1, zs = -1, zm = -1, nnodes = 0; 298 if (rank == 0) { 299 xs = grloc[r][0]; 300 xm = grloc[r][1]; 301 ys = grloc[r][2]; 302 ym = grloc[r][3]; 303 zs = grloc[r][4]; 304 zm = grloc[r][5]; 305 nnodes = xm * ym * zm; 306 } 307 maxnnodes = PetscMax(maxnnodes, nnodes); 308 PetscCall(PetscFPrintf(comm, fp, " <Piece Extent=\"%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\">\n", xs, xs + xm - 1, ys, ys + ym - 1, zs, zs + zm - 1)); 309 PetscCall(PetscFPrintf(comm, fp, " <Coordinates>\n")); 310 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"Xcoord\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, boffset)); 311 boffset += xm * sizeof(PetscScalar) + sizeof(PetscInt64); 312 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"Ycoord\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, boffset)); 313 boffset += ym * sizeof(PetscScalar) + sizeof(PetscInt64); 314 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"Zcoord\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, boffset)); 315 boffset += zm * sizeof(PetscScalar) + sizeof(PetscInt64); 316 PetscCall(PetscFPrintf(comm, fp, " </Coordinates>\n")); 317 PetscCall(PetscFPrintf(comm, fp, " <PointData Scalars=\"ScalarPointData\">\n")); 318 for (link = vtk->link; link; link = link->next) { 319 Vec X = (Vec)link->vec; 320 PetscInt bs, f; 321 DM daCurr; 322 PetscBool fieldsnamed; 323 const char *vecname = "Unnamed"; 324 325 PetscCall(VecGetDM(X, &daCurr)); 326 PetscCall(DMDAGetInfo(daCurr, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bs, NULL, NULL, NULL, NULL, NULL)); 327 maxbs = PetscMax(maxbs, bs); 328 if (((PetscObject)X)->name || link != vtk->link) PetscCall(PetscObjectGetName((PetscObject)X, &vecname)); 329 330 /* If any fields are named, add scalar fields. Otherwise, add a vector field */ 331 PetscCall(DMDAGetFieldsNamed(daCurr, &fieldsnamed)); 332 if (fieldsnamed) { 333 for (f = 0; f < bs; f++) { 334 char buf[256]; 335 const char *fieldname; 336 PetscCall(DMDAGetFieldName(daCurr, f, &fieldname)); 337 if (!fieldname) { 338 PetscCall(PetscSNPrintf(buf, sizeof(buf), "%" PetscInt_FMT, f)); 339 fieldname = buf; 340 } 341 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s.%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset)); 342 boffset += nnodes * sizeof(PetscScalar) + sizeof(PetscInt64); 343 } 344 } else { 345 PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"%" PetscInt_FMT "\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, bs, boffset)); 346 boffset += bs * nnodes * sizeof(PetscScalar) + sizeof(PetscInt64); 347 } 348 } 349 PetscCall(PetscFPrintf(comm, fp, " </PointData>\n")); 350 PetscCall(PetscFPrintf(comm, fp, " </Piece>\n")); 351 } 352 PetscCall(PetscFPrintf(comm, fp, " </RectilinearGrid>\n")); 353 PetscCall(PetscFPrintf(comm, fp, " <AppendedData encoding=\"raw\">\n")); 354 PetscCall(PetscFPrintf(comm, fp, "_")); 355 356 /* Now write the arrays. */ 357 tag = ((PetscObject)viewer)->tag; 358 PetscCall(PetscMalloc2(PetscMax(maxnnodes * maxbs, info.xm + info.ym + info.zm), &array, PetscMax(maxnnodes * maxbs, info.xm + info.ym + info.zm), &array2)); 359 for (r = 0; r < size; r++) { 360 MPI_Status status; 361 PetscInt xs = -1, xm = -1, ys = -1, ym = -1, zs = -1, zm = -1, nnodes = 0; 362 if (rank == 0) { 363 xs = grloc[r][0]; 364 xm = grloc[r][1]; 365 ys = grloc[r][2]; 366 ym = grloc[r][3]; 367 zs = grloc[r][4]; 368 zm = grloc[r][5]; 369 nnodes = xm * ym * zm; 370 } else if (r == rank) { 371 nnodes = info.xm * info.ym * info.zm; 372 } 373 { /* Write the coordinates */ 374 Vec Coords; 375 376 PetscCall(DMGetCoordinates(da, &Coords)); 377 if (Coords) { 378 const PetscScalar *coords; 379 PetscCall(VecGetArrayRead(Coords, &coords)); 380 if (rank == 0) { 381 if (r) { 382 PetscMPIInt nn; 383 PetscCallMPI(MPI_Recv(array, xm + ym + zm, MPIU_SCALAR, r, tag, comm, &status)); 384 PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &nn)); 385 PetscCheck(nn == xm + ym + zm, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Array size mismatch"); 386 } else { 387 /* x coordinates */ 388 for (j = 0, k = 0, i = 0; i < xm; i++) { 389 PetscInt Iloc = i + xm * (j + ym * k); 390 array[i] = coords[Iloc * dim + 0]; 391 } 392 /* y coordinates */ 393 for (i = 0, k = 0, j = 0; j < ym; j++) { 394 PetscInt Iloc = i + xm * (j + ym * k); 395 array[j + xm] = dim > 1 ? coords[Iloc * dim + 1] : 0; 396 } 397 /* z coordinates */ 398 for (i = 0, j = 0, k = 0; k < zm; k++) { 399 PetscInt Iloc = i + xm * (j + ym * k); 400 array[k + xm + ym] = dim > 2 ? coords[Iloc * dim + 2] : 0; 401 } 402 } 403 } else if (r == rank) { 404 xm = info.xm; 405 ym = info.ym; 406 zm = info.zm; 407 for (j = 0, k = 0, i = 0; i < xm; i++) { 408 PetscInt Iloc = i + xm * (j + ym * k); 409 array2[i] = coords[Iloc * dim + 0]; 410 } 411 for (i = 0, k = 0, j = 0; j < ym; j++) { 412 PetscInt Iloc = i + xm * (j + ym * k); 413 array2[j + xm] = dim > 1 ? coords[Iloc * dim + 1] : 0; 414 } 415 for (i = 0, j = 0, k = 0; k < zm; k++) { 416 PetscInt Iloc = i + xm * (j + ym * k); 417 array2[k + xm + ym] = dim > 2 ? coords[Iloc * dim + 2] : 0; 418 } 419 PetscCallMPI(MPI_Send((void *)array2, xm + ym + zm, MPIU_SCALAR, 0, tag, comm)); 420 } 421 PetscCall(VecRestoreArrayRead(Coords, &coords)); 422 } else { /* Fabricate some coordinates using grid index */ 423 for (i = 0; i < xm; i++) array[i] = xs + i; 424 for (j = 0; j < ym; j++) array[j + xm] = ys + j; 425 for (k = 0; k < zm; k++) array[k + xm + ym] = zs + k; 426 } 427 if (rank == 0) { 428 PetscCall(PetscViewerVTKFWrite(viewer, fp, &(array[0]), xm, MPIU_SCALAR)); 429 PetscCall(PetscViewerVTKFWrite(viewer, fp, &(array[xm]), ym, MPIU_SCALAR)); 430 PetscCall(PetscViewerVTKFWrite(viewer, fp, &(array[xm + ym]), zm, MPIU_SCALAR)); 431 } 432 } 433 434 /* Write each of the objects queued up for this file */ 435 for (link = vtk->link; link; link = link->next) { 436 Vec X = (Vec)link->vec; 437 const PetscScalar *x; 438 PetscInt bs, f; 439 DM daCurr; 440 PetscBool fieldsnamed; 441 PetscCall(VecGetDM(X, &daCurr)); 442 PetscCall(DMDAGetInfo(daCurr, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bs, NULL, NULL, NULL, NULL, NULL)); 443 444 PetscCall(VecGetArrayRead(X, &x)); 445 if (rank == 0) { 446 if (r) { 447 PetscMPIInt nn; 448 PetscCallMPI(MPI_Recv(array, nnodes * bs, MPIU_SCALAR, r, tag, comm, &status)); 449 PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &nn)); 450 PetscCheck(nn == nnodes * bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Array size mismatch receiving from rank %" PetscInt_FMT, r); 451 } else PetscCall(PetscArraycpy(array, x, nnodes * bs)); 452 /* If any fields are named, add scalar fields. Otherwise, add a vector field */ 453 PetscCall(DMDAGetFieldsNamed(daCurr, &fieldsnamed)); 454 if (fieldsnamed) { 455 for (f = 0; f < bs; f++) { 456 /* Extract and transpose the f'th field */ 457 for (k = 0; k < zm; k++) { 458 for (j = 0; j < ym; j++) { 459 for (i = 0; i < xm; i++) { 460 PetscInt Iloc = i + xm * (j + ym * k); 461 array2[Iloc] = array[Iloc * bs + f]; 462 } 463 } 464 } 465 PetscCall(PetscViewerVTKFWrite(viewer, fp, array2, nnodes, MPIU_SCALAR)); 466 } 467 } 468 PetscCall(PetscViewerVTKFWrite(viewer, fp, array, nnodes * bs, MPIU_SCALAR)); 469 470 } else if (r == rank) { 471 PetscCallMPI(MPI_Send((void *)x, nnodes * bs, MPIU_SCALAR, 0, tag, comm)); 472 } 473 PetscCall(VecRestoreArrayRead(X, &x)); 474 } 475 } 476 PetscCall(PetscFree2(array, array2)); 477 PetscCall(PetscFree(grloc)); 478 479 PetscCall(PetscFPrintf(comm, fp, "\n </AppendedData>\n")); 480 PetscCall(PetscFPrintf(comm, fp, "</VTKFile>\n")); 481 PetscCall(PetscFClose(comm, fp)); 482 PetscFunctionReturn(PETSC_SUCCESS); 483 } 484 485 /*@C 486 DMDAVTKWriteAll - Write a file containing all the fields that have been provided to the viewer 487 488 Collective 489 490 Input Parameters: 491 + odm - `DMDA` specifying the grid layout, passed as a `PetscObject` 492 - viewer - viewer of type `PETSCVIEWERVTK` 493 494 Level: developer 495 496 Notes: 497 This function is a callback used by the VTK viewer to actually write the file. 498 The reason for this odd model is that the VTK file format does not provide any way to write one field at a time. 499 Instead, metadata for the entire file needs to be available up-front before you can start writing the file. 500 501 If any fields have been named (see e.g. DMDASetFieldName()), then individual scalar 502 fields are written. Otherwise, a single multi-dof (vector) field is written. 503 504 .seealso: `DMDA`, `DM`, `PETSCVIEWERVTK` 505 @*/ 506 PetscErrorCode DMDAVTKWriteAll(PetscObject odm, PetscViewer viewer) 507 { 508 DM dm = (DM)odm; 509 PetscBool isvtk; 510 511 PetscFunctionBegin; 512 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 513 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 514 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk)); 515 PetscCheck(isvtk, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name); 516 switch (viewer->format) { 517 case PETSC_VIEWER_VTK_VTS: 518 PetscCall(DMDAVTKWriteAll_VTS(dm, viewer)); 519 break; 520 case PETSC_VIEWER_VTK_VTR: 521 PetscCall(DMDAVTKWriteAll_VTR(dm, viewer)); 522 break; 523 default: 524 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]); 525 } 526 PetscFunctionReturn(PETSC_SUCCESS); 527 } 528