1 static char help[] = "Tests adaptive refinement using DMForest, and uses HDF5.\n\n"; 2 3 #include <petscdmforest.h> 4 #include <petscdmplex.h> 5 #include <petscviewerhdf5.h> 6 7 int main (int argc, char **argv) 8 { 9 DM base, forest, plex; 10 PetscSection s; 11 PetscViewer viewer; 12 Vec g = NULL, g2 = NULL; 13 PetscReal nrm; 14 PetscBool adapt = PETSC_FALSE, userSection = PETSC_FALSE; 15 PetscInt vStart, vEnd, v, i; 16 PetscErrorCode ierr; 17 18 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 19 ierr = PetscOptionsGetBool(NULL, NULL, "-adapt", &adapt, NULL);CHKERRQ(ierr); 20 ierr = PetscOptionsGetBool(NULL, NULL, "-user_section", &userSection, NULL);CHKERRQ(ierr); 21 22 /* Create a base DMPlex mesh */ 23 ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, 2, PETSC_FALSE, NULL, NULL, NULL, NULL, PETSC_TRUE, &base);CHKERRQ(ierr); 24 ierr = DMSetFromOptions(base);CHKERRQ(ierr); 25 ierr = DMViewFromOptions(base, NULL, "-dm_view");CHKERRQ(ierr); 26 27 /* Covert Plex mesh to Forest and destroy base */ 28 ierr = DMCreate(PETSC_COMM_WORLD, &forest);CHKERRQ(ierr); 29 ierr = DMSetType(forest, DMP4EST);CHKERRQ(ierr); 30 ierr = DMForestSetBaseDM(forest, base);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 if (adapt) { 36 /* Adaptively refine the cell 0 of the mesh */ 37 for (i = 0; i < 3; ++i) { 38 DM postforest; 39 DMLabel adaptLabel = NULL; 40 41 ierr = DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel);CHKERRQ(ierr); 42 ierr = DMLabelSetValue(adaptLabel, 0, DM_ADAPT_REFINE);CHKERRQ(ierr); 43 ierr = DMForestTemplate(forest, PETSC_COMM_WORLD, &postforest);CHKERRQ(ierr); 44 ierr = DMForestSetAdaptivityLabel(postforest, adaptLabel);CHKERRQ(ierr); 45 ierr = DMLabelDestroy(&adaptLabel);CHKERRQ(ierr); 46 ierr = DMSetUp(postforest);CHKERRQ(ierr); 47 ierr = DMDestroy(&forest);CHKERRQ(ierr); 48 forest = postforest; 49 } 50 } else { 51 /* Adaptively refine all cells of the mesh */ 52 PetscInt cStart, cEnd, c; 53 54 for (i = 0; i < 3; ++i) { 55 DM postforest; 56 DMLabel adaptLabel = NULL; 57 58 ierr = DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel);CHKERRQ(ierr); 59 60 ierr = DMForestGetCellChart(forest, &cStart, &cEnd);CHKERRQ(ierr); 61 for (c = cStart; c < cEnd; ++c) { 62 ierr = DMLabelSetValue(adaptLabel, c, DM_ADAPT_REFINE);CHKERRQ(ierr); 63 } 64 65 ierr = DMForestTemplate(forest, PETSC_COMM_WORLD, &postforest);CHKERRQ(ierr); 66 ierr = DMForestSetAdaptivityLabel(postforest, adaptLabel);CHKERRQ(ierr); 67 ierr = DMLabelDestroy(&adaptLabel);CHKERRQ(ierr); 68 ierr = DMSetUp(postforest);CHKERRQ(ierr); 69 ierr = DMDestroy(&forest);CHKERRQ(ierr); 70 forest = postforest; 71 } 72 } 73 ierr = DMViewFromOptions(forest, NULL, "-dm_view");CHKERRQ(ierr); 74 75 /* Setup the section*/ 76 if (userSection) { 77 ierr = DMConvert(forest, DMPLEX, &plex);CHKERRQ(ierr); 78 ierr = DMViewFromOptions(plex, NULL, "-dm_view");CHKERRQ(ierr); 79 ierr = DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd);CHKERRQ(ierr); 80 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "Vertices [%D, %D)\n", vStart, vEnd);CHKERRQ(ierr); 81 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL);CHKERRQ(ierr); 82 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) forest), &s);CHKERRQ(ierr); 83 ierr = PetscSectionSetNumFields(s, 1);CHKERRQ(ierr); 84 ierr = PetscSectionSetChart(s, vStart, vEnd);CHKERRQ(ierr); 85 for (v = vStart; v < vEnd; ++v) { 86 ierr = PetscSectionSetDof(s, v, 1);CHKERRQ(ierr); 87 ierr = PetscSectionSetFieldDof(s, v, 0, 1);CHKERRQ(ierr); 88 } 89 ierr = PetscSectionSetUp(s);CHKERRQ(ierr); 90 ierr = DMSetLocalSection(forest, s);CHKERRQ(ierr); 91 ierr = PetscObjectViewFromOptions((PetscObject) s, NULL, "-my_section_view");CHKERRQ(ierr); 92 ierr = PetscSectionDestroy(&s);CHKERRQ(ierr); 93 ierr = DMDestroy(&plex);CHKERRQ(ierr); 94 } else { 95 PetscFE fe; 96 PetscInt dim; 97 98 ierr = DMGetDimension(forest, &dim);CHKERRQ(ierr); 99 ierr = PetscFECreateLagrange(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, 1, PETSC_DETERMINE, &fe);CHKERRQ(ierr); 100 ierr = DMAddField(forest, NULL, (PetscObject) fe);CHKERRQ(ierr); 101 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 102 ierr = DMCreateDS(forest);CHKERRQ(ierr); 103 } 104 105 /* Create the global vector*/ 106 ierr = DMCreateGlobalVector(forest, &g);CHKERRQ(ierr); 107 ierr = PetscObjectSetName((PetscObject) g, "g");CHKERRQ(ierr); 108 ierr = VecSet(g, 1.0);CHKERRQ(ierr); 109 110 /* Test global to local*/ 111 Vec l; 112 ierr = DMCreateLocalVector(forest, &l);CHKERRQ(ierr); 113 ierr = VecZeroEntries(l);CHKERRQ(ierr); 114 ierr = DMGlobalToLocal(forest, g, INSERT_VALUES, l);CHKERRQ(ierr); 115 ierr = VecZeroEntries(g);CHKERRQ(ierr); 116 ierr = DMLocalToGlobal(forest, l, INSERT_VALUES, g);CHKERRQ(ierr); 117 ierr = VecDestroy(&l);CHKERRQ(ierr); 118 119 /* Save a vector*/ 120 ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, "forestHDF.h5", FILE_MODE_WRITE, &viewer); 121 ierr = VecView(g, viewer);CHKERRQ(ierr); 122 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 123 124 /* Load another vector to load into*/ 125 ierr = DMCreateGlobalVector(forest, &g2);CHKERRQ(ierr); 126 ierr = PetscObjectSetName((PetscObject) g2, "g");CHKERRQ(ierr); 127 ierr = VecZeroEntries(g2);CHKERRQ(ierr); 128 129 /* Load a vector*/ 130 ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, "forestHDF.h5", FILE_MODE_READ, &viewer); 131 ierr = VecLoad(g2, viewer);CHKERRQ(ierr); 132 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 133 134 /* Check if the data is the same*/ 135 ierr = VecAXPY(g2, -1.0, g);CHKERRQ(ierr); 136 ierr = VecNorm(g2, NORM_INFINITY, &nrm);CHKERRQ(ierr); 137 if (PetscAbsReal(nrm) > PETSC_SMALL) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Invalid difference norm %g", (double) nrm); 138 139 ierr = VecDestroy(&g);CHKERRQ(ierr); 140 ierr = VecDestroy(&g2);CHKERRQ(ierr); 141 ierr = DMDestroy(&forest);CHKERRQ(ierr); 142 ierr = PetscFinalize(); 143 return ierr; 144 } 145 146 /*TEST 147 148 build: 149 requires: hdf5 p4est 150 151 test: 152 suffix: 0 153 nsize: {{1 2 5}} 154 args: -adapt -dm_plex_box_faces 2,2 155 156 TEST*/ 157