1 static char help[] = "Test HDF5 input and output.\n\ 2 This exposed a bug with sharing discretizations.\n\n\n"; 3 4 #include <petscdmforest.h> 5 #include <petscdmplex.h> 6 #include <petscviewerhdf5.h> 7 8 int main(int argc, char **argv) { 9 DM base, forest, plex; 10 Vec g, g2; 11 PetscSection s; 12 PetscViewer viewer; 13 PetscReal diff; 14 PetscInt min_refine = 2, overlap = 0; 15 PetscInt vStart, vEnd, v; 16 17 PetscFunctionBeginUser; 18 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 19 20 PetscCall(DMCreate(PETSC_COMM_WORLD, &base)); 21 PetscCall(DMSetType(base, DMPLEX)); 22 PetscCall(DMSetFromOptions(base)); 23 24 PetscCall(DMCreate(PETSC_COMM_WORLD, &forest)); 25 PetscCall(DMSetType(forest, DMP4EST)); 26 PetscCall(DMForestSetBaseDM(forest, base)); 27 PetscCall(DMForestSetInitialRefinement(forest, min_refine)); 28 PetscCall(DMForestSetPartitionOverlap(forest, overlap)); 29 PetscCall(DMSetUp(forest)); 30 PetscCall(DMDestroy(&base)); 31 PetscCall(DMViewFromOptions(forest, NULL, "-dm_view")); 32 33 PetscCall(DMConvert(forest, DMPLEX, &plex)); 34 PetscCall(DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd)); 35 PetscCall(DMDestroy(&plex)); 36 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)forest), &s)); 37 PetscCall(PetscSectionSetChart(s, vStart, vEnd)); 38 for (v = vStart; v < vEnd; ++v) PetscCall(PetscSectionSetDof(s, v, 1)); 39 PetscCall(PetscSectionSetUp(s)); 40 PetscCall(DMSetLocalSection(forest, s)); 41 PetscCall(PetscSectionDestroy(&s)); 42 43 PetscCall(DMCreateGlobalVector(forest, &g)); 44 PetscCall(PetscObjectSetName((PetscObject)g, "g")); 45 PetscCall(VecSet(g, 1.0)); 46 PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "forest.h5", FILE_MODE_WRITE, &viewer)); 47 PetscCall(VecView(g, viewer)); 48 PetscCall(PetscViewerDestroy(&viewer)); 49 50 PetscCall(DMCreateGlobalVector(forest, &g2)); 51 PetscCall(PetscObjectSetName((PetscObject)g2, "g")); 52 PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "forest.h5", FILE_MODE_READ, &viewer)); 53 PetscCall(VecLoad(g2, viewer)); 54 PetscCall(PetscViewerDestroy(&viewer)); 55 56 /* Check if the data is the same*/ 57 PetscCall(VecAXPY(g2, -1.0, g)); 58 PetscCall(VecNorm(g2, NORM_INFINITY, &diff)); 59 if (diff > PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Check failed: %g\n", (double)diff)); 60 61 PetscCall(VecDestroy(&g)); 62 PetscCall(VecDestroy(&g2)); 63 PetscCall(DMDestroy(&forest)); 64 PetscCall(PetscFinalize()); 65 return 0; 66 } 67 68 /*TEST 69 70 build: 71 requires: hdf5 p4est 72 73 test: 74 args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3 75 76 TEST*/ 77