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 PetscCheck(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 PetscCheck(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 PetscCheck(!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 314 CHKERRQ(PetscInitialize(&argc, &argv, NULL, help)); 315 CHKERRQ(ProcessOptions(PETSC_COMM_WORLD, &user)); 316 CHKERRQ(CreateMesh(&user, &dm)); 317 if (user.compare_boundary) { 318 CHKERRQ(DMAddBoundaryLabel_GetCoordinateRepresentation(dm, &allCoords)); 319 } 320 CHKERRQ(SaveMesh(&user, dm)); 321 CHKERRQ(LoadMesh(&user, &dmnew)); 322 CHKERRQ(CompareMeshes(&user, dm, dmnew)); 323 if (user.compare_boundary) { 324 CHKERRQ(DMGetBoundaryLabel_CompareWithCoordinateRepresentation(&user, dmnew, allCoords)); 325 } 326 CHKERRQ(DMDestroy(&dm)); 327 CHKERRQ(DMDestroy(&dmnew)); 328 CHKERRQ(VecDestroy(&allCoords)); 329 CHKERRQ(PetscFinalize()); 330 return 0; 331 } 332 333 //TODO we can -compare once the new parallel topology format is in place 334 /*TEST 335 build: 336 requires: hdf5 337 338 # load old format, save in new format, reload, distribute 339 testset: 340 suffix: 1 341 requires: !complex datafilespath 342 args: -dm_plex_name plex 343 args: -dm_plex_check_all -dm_plex_view_hdf5_storage_version 2.0.0 344 args: -dm_plex_interpolate 345 args: -load_dm_plex_check_all 346 args: -use_low_level_functions {{0 1}} -compare_boundary 347 args: -outfile ex56_1.h5 348 nsize: {{1 3}} 349 test: 350 suffix: a 351 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/annulus-20.h5 352 test: 353 suffix: b 354 TODO: broken 355 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/barycentricallyrefinedcube.h5 356 test: 357 suffix: c 358 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/blockcylinder-50.h5 359 test: 360 suffix: d 361 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/cube-hexahedra-refined.h5 362 test: 363 suffix: e 364 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/hybrid_hexwedge.h5 365 test: 366 suffix: f 367 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/square.h5 368 369 # load old format, save in new format, reload topology, distribute, load geometry and labels 370 testset: 371 suffix: 2 372 requires: !complex datafilespath 373 args: -dm_plex_name plex 374 args: -dm_plex_check_all -dm_plex_view_hdf5_storage_version 2.0.0 375 args: -dm_plex_interpolate 376 args: -load_dm_plex_check_all 377 args: -use_low_level_functions -load_dm_distribute 0 -distribute_after_topo_load -compare_boundary 378 args: -outfile ex56_2.h5 379 nsize: 3 380 test: 381 suffix: a 382 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/annulus-20.h5 383 test: 384 suffix: b 385 TODO: broken 386 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/barycentricallyrefinedcube.h5 387 test: 388 suffix: c 389 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/blockcylinder-50.h5 390 test: 391 suffix: d 392 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/cube-hexahedra-refined.h5 393 test: 394 suffix: e 395 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/hybrid_hexwedge.h5 396 test: 397 suffix: f 398 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/square.h5 399 400 # load old format, save in new format, reload topology, distribute, load geometry and labels 401 testset: 402 suffix: 3 403 requires: !complex datafilespath 404 args: -dm_plex_name plex 405 args: -dm_plex_view_hdf5_storage_version 2.0.0 406 args: -dm_plex_interpolate -load_dm_distribute 0 407 args: -use_low_level_functions -compare_pre_post 408 args: -outfile ex56_3.h5 409 nsize: 3 410 test: 411 suffix: a 412 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/annulus-20.h5 413 test: 414 suffix: b 415 TODO: broken 416 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/barycentricallyrefinedcube.h5 417 test: 418 suffix: c 419 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/blockcylinder-50.h5 420 test: 421 suffix: d 422 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/cube-hexahedra-refined.h5 423 test: 424 suffix: e 425 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/hybrid_hexwedge.h5 426 test: 427 suffix: f 428 args: -dm_plex_filename ${DATAFILESPATH}/meshes/hdf5-petsc/petsc-v3.16.0/v1.0.0/square.h5 429 TEST*/ 430