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