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