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 ierr = DMCreate(PETSC_COMM_WORLD, &base);CHKERRQ(ierr); 23 ierr = DMSetType(base, DMPLEX);CHKERRQ(ierr); 24 ierr = DMSetFromOptions(base);CHKERRQ(ierr); 25 26 ierr = DMCreate(PETSC_COMM_WORLD, &forest);CHKERRQ(ierr); 27 ierr = DMSetType(forest, DMP4EST);CHKERRQ(ierr); 28 ierr = DMForestSetBaseDM(forest, base);CHKERRQ(ierr); 29 ierr = DMForestSetInitialRefinement(forest, min_refine);CHKERRQ(ierr); 30 ierr = DMForestSetPartitionOverlap(forest, overlap);CHKERRQ(ierr); 31 ierr = DMSetUp(forest);CHKERRQ(ierr); 32 ierr = DMDestroy(&base);CHKERRQ(ierr); 33 ierr = DMViewFromOptions(forest, NULL, "-dm_view");CHKERRQ(ierr); 34 35 ierr = DMConvert(forest, DMPLEX, &plex);CHKERRQ(ierr); 36 ierr = DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd);CHKERRQ(ierr); 37 ierr = DMDestroy(&plex);CHKERRQ(ierr); 38 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) forest), &s);CHKERRQ(ierr); 39 ierr = PetscSectionSetChart(s, vStart, vEnd);CHKERRQ(ierr); 40 for (v = vStart; v < vEnd; ++v) {ierr = PetscSectionSetDof(s, v, 1);CHKERRQ(ierr);} 41 ierr = PetscSectionSetUp(s);CHKERRQ(ierr); 42 ierr = DMSetLocalSection(forest, s);CHKERRQ(ierr); 43 ierr = PetscSectionDestroy(&s);CHKERRQ(ierr); 44 45 ierr = DMCreateGlobalVector(forest, &g);CHKERRQ(ierr); 46 ierr = PetscObjectSetName((PetscObject) g, "g");CHKERRQ(ierr); 47 ierr = VecSet(g, 1.0);CHKERRQ(ierr); 48 ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, "forest.h5", FILE_MODE_WRITE, &viewer); 49 ierr = VecView(g, viewer);CHKERRQ(ierr); 50 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 51 52 ierr = DMCreateGlobalVector(forest, &g2);CHKERRQ(ierr); 53 ierr = PetscObjectSetName((PetscObject) g2, "g");CHKERRQ(ierr); 54 ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, "forest.h5", FILE_MODE_READ, &viewer); 55 ierr = VecLoad(g2, viewer);CHKERRQ(ierr); 56 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 57 58 /* Check if the data is the same*/ 59 ierr = VecAXPY(g2, -1.0, g);CHKERRQ(ierr); 60 ierr = VecNorm(g2, NORM_INFINITY, &diff);CHKERRQ(ierr); 61 if (diff > PETSC_MACHINE_EPSILON) {ierr = PetscPrintf(PETSC_COMM_WORLD, "Check failed: %g\n", (double) diff);CHKERRQ(ierr);} 62 63 ierr = VecDestroy(&g);CHKERRQ(ierr); 64 ierr = VecDestroy(&g2);CHKERRQ(ierr); 65 ierr = DMDestroy(&forest);CHKERRQ(ierr); 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