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