1 #include <petscdmplex.h> 2 #include <petscviewerhdf5.h> 3 #include <petscsf.h> 4 5 static const char help[] = "Load and save the mesh to the native HDF5 format\n\n"; 6 static const char EX[] = "ex56.c"; 7 static const char LABEL_NAME[] = "BoundaryVertices"; 8 static const PetscInt LABEL_VALUE = 12345; 9 typedef struct { 10 MPI_Comm comm; 11 const char *meshname; /* Mesh name */ 12 PetscBool compare; /* Compare the meshes using DMPlexEqual() and DMCompareLabels() */ 13 PetscBool compare_labels; /* Compare labels in the meshes using DMCompareLabels() */ 14 PetscBool compare_boundary; /* Check label I/O via boundary vertex coordinates */ 15 PetscBool compare_pre_post; /* Compare labels loaded before distribution with those loaded after distribution */ 16 char outfile[PETSC_MAX_PATH_LEN]; /* Output file */ 17 PetscBool use_low_level_functions; /* Use low level functions for viewing and loading */ 18 //TODO This is meant as temporary option; can be removed once we have full parallel loading in place 19 PetscBool distribute_after_topo_load; /* Distribute topology right after DMPlexTopologyLoad(), if use_low_level_functions=true */ 20 PetscBool verbose; 21 } AppCtx; 22 23 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 24 { 25 PetscErrorCode ierr; 26 27 PetscFunctionBeginUser; 28 options->comm = comm; 29 options->compare = PETSC_FALSE; 30 options->compare_labels = PETSC_FALSE; 31 options->compare_boundary = PETSC_FALSE; 32 options->compare_pre_post = PETSC_FALSE; 33 options->outfile[0] = '\0'; 34 options->use_low_level_functions = PETSC_FALSE; 35 options->distribute_after_topo_load = PETSC_FALSE; 36 options->verbose = PETSC_FALSE; 37 38 ierr = PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");CHKERRQ(ierr); 39 CHKERRQ(PetscOptionsBool("-compare", "Compare the meshes using DMPlexEqual() and DMCompareLabels()", EX, options->compare, &options->compare, NULL)); 40 CHKERRQ(PetscOptionsBool("-compare_labels", "Compare labels in the meshes using DMCompareLabels()", "ex55.c", options->compare_labels, &options->compare_labels, NULL)); 41 CHKERRQ(PetscOptionsBool("-compare_boundary", "Check label I/O via boundary vertex coordinates", "ex55.c", options->compare_boundary, &options->compare_boundary, NULL)); 42 CHKERRQ(PetscOptionsBool("-compare_pre_post", "Compare labels loaded before distribution with those loaded after distribution", "ex55.c", options->compare_pre_post, &options->compare_pre_post, NULL)); 43 CHKERRQ(PetscOptionsString("-outfile", "Output mesh file", EX, options->outfile, options->outfile, sizeof(options->outfile), NULL)); 44 CHKERRQ(PetscOptionsBool("-use_low_level_functions", "Use low level functions for viewing and loading", EX, options->use_low_level_functions, &options->use_low_level_functions, NULL)); 45 CHKERRQ(PetscOptionsBool("-distribute_after_topo_load", "Distribute topology right after DMPlexTopologyLoad(), if use_low_level_functions=true", EX, options->distribute_after_topo_load, &options->distribute_after_topo_load, NULL)); 46 CHKERRQ(PetscOptionsBool("-verbose", "Verbose printing", EX, options->verbose, &options->verbose, NULL)); 47 ierr = PetscOptionsEnd();CHKERRQ(ierr); 48 PetscFunctionReturn(0); 49 }; 50 51 static PetscErrorCode CreateMesh(AppCtx *options, DM *newdm) 52 { 53 DM dm; 54 55 PetscFunctionBeginUser; 56 CHKERRQ(DMCreate(options->comm, &dm)); 57 CHKERRQ(DMSetType(dm, DMPLEX)); 58 CHKERRQ(DMSetFromOptions(dm)); 59 CHKERRQ(PetscObjectGetName((PetscObject)dm, &options->meshname)); 60 CHKERRQ(DMViewFromOptions(dm, NULL, "-dm_view")); 61 *newdm = dm; 62 PetscFunctionReturn(0); 63 } 64 65 static PetscErrorCode SaveMesh(AppCtx *options, DM dm) 66 { 67 PetscViewer v; 68 69 PetscFunctionBeginUser; 70 CHKERRQ(PetscViewerHDF5Open(PetscObjectComm((PetscObject) dm), options->outfile, FILE_MODE_WRITE, &v)); 71 if (options->use_low_level_functions) { 72 CHKERRQ(DMPlexTopologyView(dm, v)); 73 CHKERRQ(DMPlexCoordinatesView(dm, v)); 74 CHKERRQ(DMPlexLabelsView(dm, v)); 75 } else { 76 CHKERRQ(DMView(dm, v)); 77 } 78 CHKERRQ(PetscViewerDestroy(&v)); 79 PetscFunctionReturn(0); 80 } 81 82 typedef enum {NONE=0, PRE_DIST=1, POST_DIST=2} AuxObjLoadMode; 83 84 static PetscErrorCode LoadMeshLowLevel(AppCtx *options, PetscViewer v, PetscBool explicitDistribute, AuxObjLoadMode mode, DM *newdm) 85 { 86 DM dm; 87 PetscSF sfXC; 88 89 PetscFunctionBeginUser; 90 CHKERRQ(DMCreate(options->comm, &dm)); 91 CHKERRQ(DMSetType(dm, DMPLEX)); 92 CHKERRQ(PetscObjectSetName((PetscObject) dm, options->meshname)); 93 CHKERRQ(DMPlexTopologyLoad(dm, v, &sfXC)); 94 if (mode == PRE_DIST) { 95 CHKERRQ(DMPlexCoordinatesLoad(dm, v, sfXC)); 96 CHKERRQ(DMPlexLabelsLoad(dm, v, sfXC)); 97 } 98 if (explicitDistribute) { 99 DM dmdist; 100 PetscSF sfXB = sfXC, sfBC; 101 102 CHKERRQ(DMPlexDistribute(dm, 0, &sfBC, &dmdist)); 103 if (dmdist) { 104 const char *name; 105 106 CHKERRQ(PetscObjectGetName((PetscObject) dm, &name)); 107 CHKERRQ(PetscObjectSetName((PetscObject) dmdist, name)); 108 CHKERRQ(PetscSFCompose(sfXB, sfBC, &sfXC)); 109 CHKERRQ(PetscSFDestroy(&sfXB)); 110 CHKERRQ(PetscSFDestroy(&sfBC)); 111 CHKERRQ(DMDestroy(&dm)); 112 dm = dmdist; 113 } 114 } 115 if (mode == POST_DIST) { 116 CHKERRQ(DMPlexCoordinatesLoad(dm, v, sfXC)); 117 CHKERRQ(DMPlexLabelsLoad(dm, v, sfXC)); 118 } 119 CHKERRQ(PetscSFDestroy(&sfXC)); 120 *newdm = dm; 121 PetscFunctionReturn(0); 122 } 123 124 static PetscErrorCode LoadMesh(AppCtx *options, DM *dmnew) 125 { 126 DM dm; 127 PetscViewer v; 128 129 PetscFunctionBeginUser; 130 CHKERRQ(PetscViewerHDF5Open(options->comm, options->outfile, FILE_MODE_READ, &v)); 131 if (options->use_low_level_functions) { 132 if (options->compare_pre_post) { 133 DM dm0; 134 135 CHKERRQ(LoadMeshLowLevel(options, v, PETSC_TRUE, PRE_DIST, &dm0)); 136 CHKERRQ(LoadMeshLowLevel(options, v, PETSC_TRUE, POST_DIST, &dm)); 137 CHKERRQ(DMCompareLabels(dm0, dm, NULL, NULL)); 138 CHKERRQ(DMDestroy(&dm0)); 139 } else { 140 CHKERRQ(LoadMeshLowLevel(options, v, options->distribute_after_topo_load, POST_DIST, &dm)); 141 } 142 } else { 143 CHKERRQ(DMCreate(options->comm, &dm)); 144 CHKERRQ(DMSetType(dm, DMPLEX)); 145 CHKERRQ(PetscObjectSetName((PetscObject) dm, options->meshname)); 146 CHKERRQ(DMLoad(dm, v)); 147 } 148 CHKERRQ(PetscViewerDestroy(&v)); 149 150 CHKERRQ(DMSetOptionsPrefix(dm, "load_")); 151 CHKERRQ(DMSetFromOptions(dm)); 152 CHKERRQ(DMViewFromOptions(dm, NULL, "-dm_view")); 153 *dmnew = dm; 154 PetscFunctionReturn(0); 155 } 156 157 static PetscErrorCode CompareMeshes(AppCtx *options, DM dm0, DM dm1) 158 { 159 PetscBool flg; 160 161 PetscFunctionBeginUser; 162 if (options->compare) { 163 CHKERRQ(DMPlexEqual(dm0, dm1, &flg)); 164 PetscCheckFalse(!flg,options->comm, PETSC_ERR_ARG_INCOMP, "DMs are not equal"); 165 CHKERRQ(PetscPrintf(options->comm,"DMs equal\n")); 166 } 167 if (options->compare_labels) { 168 CHKERRQ(DMCompareLabels(dm0, dm1, NULL, NULL)); 169 CHKERRQ(PetscPrintf(options->comm,"DMLabels equal\n")); 170 } 171 PetscFunctionReturn(0); 172 } 173 174 static PetscErrorCode MarkBoundaryVertices(DM dm, PetscInt value, DMLabel *label) 175 { 176 DMLabel l; 177 IS points; 178 const PetscInt *idx; 179 PetscInt i, n; 180 181 PetscFunctionBeginUser; 182 CHKERRQ(DMLabelCreate(PetscObjectComm((PetscObject)dm), LABEL_NAME, &l)); 183 CHKERRQ(DMPlexMarkBoundaryFaces(dm, value, l)); 184 CHKERRQ(DMPlexLabelComplete(dm, l)); 185 CHKERRQ(DMLabelGetStratumIS(l, value, &points)); 186 187 CHKERRQ(ISGetLocalSize(points, &n)); 188 CHKERRQ(ISGetIndices(points, &idx)); 189 for (i=0; i<n; i++) { 190 const PetscInt p = idx[i]; 191 PetscInt d; 192 193 CHKERRQ(DMPlexGetPointDepth(dm, p, &d)); 194 if (d != 0) { 195 CHKERRQ(DMLabelClearValue(l, p, value)); 196 } 197 } 198 CHKERRQ(ISRestoreIndices(points, &idx)); 199 CHKERRQ(ISDestroy(&points)); 200 *label = l; 201 PetscFunctionReturn(0); 202 } 203 204 static PetscErrorCode VertexCoordinatesToAll(DM dm, IS vertices, Vec *allCoords) 205 { 206 Vec coords, allCoords_; 207 VecScatter sc; 208 MPI_Comm comm; 209 210 PetscFunctionBeginUser; 211 CHKERRQ(PetscObjectGetComm((PetscObject)dm, &comm)); 212 CHKERRQ(DMGetCoordinatesLocalSetUp(dm)); 213 if (vertices) { 214 CHKERRQ(DMGetCoordinatesLocalTuple(dm, vertices, NULL, &coords)); 215 } else { 216 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF, 0, &coords)); 217 } 218 { 219 PetscInt n; 220 Vec mpivec; 221 222 CHKERRQ(VecGetLocalSize(coords, &n)); 223 CHKERRQ(VecCreateMPI(comm, n, PETSC_DECIDE, &mpivec)); 224 CHKERRQ(VecCopy(coords, mpivec)); 225 CHKERRQ(VecDestroy(&coords)); 226 coords = mpivec; 227 } 228 229 CHKERRQ(VecScatterCreateToAll(coords, &sc, &allCoords_)); 230 CHKERRQ(VecScatterBegin(sc,coords,allCoords_,INSERT_VALUES,SCATTER_FORWARD)); 231 CHKERRQ(VecScatterEnd(sc,coords,allCoords_,INSERT_VALUES,SCATTER_FORWARD)); 232 CHKERRQ(VecScatterDestroy(&sc)); 233 CHKERRQ(VecDestroy(&coords)); 234 *allCoords = allCoords_; 235 PetscFunctionReturn(0); 236 } 237 238 /* Compute boundary label, remember boundary vertices using coordinates, save and load label, check it is defined on the original boundary vertices */ 239 static PetscErrorCode DMAddBoundaryLabel_GetCoordinateRepresentation(DM dm, Vec *allCoords) 240 { 241 DMLabel label; 242 IS vertices; 243 244 PetscFunctionBeginUser; 245 CHKERRQ(MarkBoundaryVertices(dm, LABEL_VALUE, &label)); 246 CHKERRQ(DMLabelGetStratumIS(label, LABEL_VALUE, &vertices)); 247 CHKERRQ(VertexCoordinatesToAll(dm, vertices, allCoords)); 248 CHKERRQ(DMAddLabel(dm, label)); 249 CHKERRQ(DMLabelDestroy(&label)); 250 CHKERRQ(ISDestroy(&vertices)); 251 PetscFunctionReturn(0); 252 } 253 254 static PetscErrorCode DMGetBoundaryLabel_CompareWithCoordinateRepresentation(AppCtx *user, DM dm, Vec allCoords) 255 { 256 DMLabel label; 257 IS pointsIS; 258 const PetscInt *points; 259 PetscInt i, n; 260 PetscBool fail = PETSC_FALSE; 261 MPI_Comm comm; 262 PetscMPIInt rank; 263 264 PetscFunctionBeginUser; 265 CHKERRQ(PetscObjectGetComm((PetscObject)dm, &comm)); 266 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 267 CHKERRQ(DMGetLabel(dm, LABEL_NAME, &label)); 268 PetscCheckFalse(!label,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Label \"%s\" was not loaded", LABEL_NAME); 269 { 270 PetscInt pStart, pEnd; 271 272 CHKERRQ(DMPlexGetChart(dm, &pStart, &pEnd)); 273 CHKERRQ(DMLabelCreateIndex(label, pStart, pEnd)); 274 } 275 CHKERRQ(DMPlexFindVertices(dm, allCoords, 0.0, &pointsIS)); 276 CHKERRQ(ISGetIndices(pointsIS, &points)); 277 CHKERRQ(ISGetLocalSize(pointsIS, &n)); 278 if (user->verbose) CHKERRQ(DMLabelView(label, PETSC_VIEWER_STDOUT_(comm))); 279 for (i=0; i<n; i++) { 280 const PetscInt p = points[i]; 281 PetscBool has; 282 PetscInt v; 283 284 if (p < 0) continue; 285 CHKERRQ(DMLabelHasPoint(label, p, &has)); 286 if (!has) { 287 CHKERRQ(PetscSynchronizedFPrintf(comm, PETSC_STDERR, "[%d] Label does not have point %D\n", rank, p)); 288 fail = PETSC_TRUE; 289 continue; 290 } 291 CHKERRQ(DMLabelGetValue(label, p, &v)); 292 if (v != LABEL_VALUE) { 293 CHKERRQ(PetscSynchronizedFPrintf(comm, PETSC_STDERR, "Point %D has bad value %D", p, v)); 294 fail = PETSC_TRUE; 295 continue; 296 } 297 if (user->verbose) CHKERRQ(PetscSynchronizedPrintf(comm, "[%d] OK point %D\n", rank, p)); 298 } 299 CHKERRQ(PetscSynchronizedFlush(comm, PETSC_STDOUT)); 300 CHKERRQ(PetscSynchronizedFlush(comm, PETSC_STDERR)); 301 CHKERRQ(ISRestoreIndices(pointsIS, &points)); 302 CHKERRQ(ISDestroy(&pointsIS)); 303 CHKERRMPI(MPI_Allreduce(MPI_IN_PLACE, &fail, 1, MPIU_BOOL, MPI_LOR, comm)); 304 PetscCheckFalse(fail,comm, PETSC_ERR_PLIB, "Label \"%s\" was not loaded correctly - see details above", LABEL_NAME); 305 PetscFunctionReturn(0); 306 } 307 308 int main(int argc, char **argv) 309 { 310 DM dm, dmnew; 311 AppCtx user; 312 Vec allCoords = NULL; 313 PetscErrorCode ierr; 314 315 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 316 CHKERRQ(ProcessOptions(PETSC_COMM_WORLD, &user)); 317 CHKERRQ(CreateMesh(&user, &dm)); 318 if (user.compare_boundary) { 319 CHKERRQ(DMAddBoundaryLabel_GetCoordinateRepresentation(dm, &allCoords)); 320 } 321 CHKERRQ(SaveMesh(&user, dm)); 322 CHKERRQ(LoadMesh(&user, &dmnew)); 323 CHKERRQ(CompareMeshes(&user, dm, dmnew)); 324 if (user.compare_boundary) { 325 CHKERRQ(DMGetBoundaryLabel_CompareWithCoordinateRepresentation(&user, dmnew, allCoords)); 326 } 327 CHKERRQ(DMDestroy(&dm)); 328 CHKERRQ(DMDestroy(&dmnew)); 329 CHKERRQ(VecDestroy(&allCoords)); 330 ierr = PetscFinalize(); 331 return ierr; 332 } 333 334 //TODO we can -compare once the new parallel topology format is in place 335 /*TEST 336 build: 337 requires: hdf5 338 339 # load old format, save in new format, reload, distribute 340 testset: 341 suffix: 1 342 requires: !complex datafilespath 343 args: -dm_plex_name plex 344 args: -dm_plex_check_all -dm_plex_view_hdf5_storage_version 2.0.0 345 args: -dm_plex_interpolate 346 args: -load_dm_plex_check_all 347 args: -use_low_level_functions {{0 1}} -compare_boundary 348 args: -outfile ex56_1.h5 349 nsize: {{1 3}} 350 test: 351 suffix: a 352 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/annulus-20.h5 353 test: 354 suffix: b 355 TODO: broken 356 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/barycentricallyrefinedcube.h5 357 test: 358 suffix: c 359 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/blockcylinder-50.h5 360 test: 361 suffix: d 362 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/cube-hexahedra-refined.h5 363 test: 364 suffix: e 365 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/hybrid_hexwedge.h5 366 test: 367 suffix: f 368 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/square.h5 369 370 # load old format, save in new format, reload topology, distribute, load geometry and labels 371 testset: 372 suffix: 2 373 requires: !complex datafilespath 374 args: -dm_plex_name plex 375 args: -dm_plex_check_all -dm_plex_view_hdf5_storage_version 2.0.0 376 args: -dm_plex_interpolate 377 args: -load_dm_plex_check_all 378 args: -use_low_level_functions -load_dm_distribute 0 -distribute_after_topo_load -compare_boundary 379 args: -outfile ex56_2.h5 380 nsize: 3 381 test: 382 suffix: a 383 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/annulus-20.h5 384 test: 385 suffix: b 386 TODO: broken 387 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/barycentricallyrefinedcube.h5 388 test: 389 suffix: c 390 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/blockcylinder-50.h5 391 test: 392 suffix: d 393 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/cube-hexahedra-refined.h5 394 test: 395 suffix: e 396 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/hybrid_hexwedge.h5 397 test: 398 suffix: f 399 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/square.h5 400 401 # load old format, save in new format, reload topology, distribute, load geometry and labels 402 testset: 403 suffix: 3 404 requires: !complex datafilespath 405 args: -dm_plex_name plex 406 args: -dm_plex_view_hdf5_storage_version 2.0.0 407 args: -dm_plex_interpolate -load_dm_distribute 0 408 args: -use_low_level_functions -compare_pre_post 409 args: -outfile ex56_3.h5 410 nsize: 3 411 test: 412 suffix: a 413 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/annulus-20.h5 414 test: 415 suffix: b 416 TODO: broken 417 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/barycentricallyrefinedcube.h5 418 test: 419 suffix: c 420 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/blockcylinder-50.h5 421 test: 422 suffix: d 423 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/cube-hexahedra-refined.h5 424 test: 425 suffix: e 426 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/hybrid_hexwedge.h5 427 test: 428 suffix: f 429 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/square.h5 430 TEST*/ 431