1 static char help[] = "TDycore Mesh Examples\n\n"; 2 3 #include <petscdmplex.h> 4 5 typedef struct { 6 char filename[PETSC_MAX_PATH_LEN]; /* Import mesh from file */ 7 PetscBool adapt; /* Flag for adaptation of the surface mesh */ 8 PetscBool extrude; /* Flag for extrusion of the suraace mesh */ 9 } AppCtx; 10 11 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 12 { 13 PetscErrorCode ierr; 14 15 PetscFunctionBeginUser; 16 options->filename[0] = '\0'; 17 options->adapt = PETSC_FALSE; 18 options->extrude = PETSC_TRUE; 19 20 ierr = PetscOptionsBegin(comm, "", "Meshing Interpolation Test Options", "DMPLEX");CHKERRQ(ierr); 21 ierr = PetscOptionsString("-filename", "The mesh file", "ex10.c", options->filename, options->filename, sizeof(options->filename), NULL);CHKERRQ(ierr); 22 ierr = PetscOptionsBool("-adapt", "Flag for adaptation of the surface mesh", "ex10.c", options->adapt, &options->adapt, NULL);CHKERRQ(ierr); 23 ierr = PetscOptionsBool("-extrude", "Flag for extrusion of the surface mesh", "ex10.c", options->extrude, &options->extrude, NULL);CHKERRQ(ierr); 24 ierr = PetscOptionsEnd(); 25 PetscFunctionReturn(0); 26 } 27 28 static PetscErrorCode CreateDomainLabel(DM dm) 29 { 30 DMLabel label; 31 PetscInt cStart, cEnd, c; 32 PetscErrorCode ierr; 33 34 PetscFunctionBeginUser; 35 ierr = DMCreateLabel(dm, "Cell Sets");CHKERRQ(ierr); 36 ierr = DMGetLabel(dm, "Cell Sets", &label);CHKERRQ(ierr); 37 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 38 for (c = cStart; c < cEnd; ++c) { 39 PetscReal centroid[3], volume, x, y; 40 41 ierr = DMPlexComputeCellGeometryFVM(dm, c, &volume, centroid, NULL);CHKERRQ(ierr); 42 x = centroid[0]; y = centroid[1]; 43 /* Headwaters are (0.0,0.25)--(0.1,0.75) */ 44 if ((x >= 0.0 && x < 0.1) && (y >= 0.25 && y <= 0.75)) {ierr = DMLabelSetValue(label, c, 1);CHKERRQ(ierr);continue;} 45 /* River channel is (0.1,0.45)--(1.0,0.55) */ 46 if ((x >= 0.1 && x <= 1.0) && (y >= 0.45 && y <= 0.55)) {ierr = DMLabelSetValue(label, c, 2);CHKERRQ(ierr);continue;} 47 } 48 PetscFunctionReturn(0); 49 } 50 51 static PetscErrorCode AdaptMesh(DM *dm, AppCtx *ctx) 52 { 53 DM dmCur = *dm; 54 DMLabel label; 55 IS valueIS, vIS; 56 PetscBool hasLabel; 57 const PetscInt *values; 58 PetscReal *volConst; /* Volume constraints for each label value */ 59 PetscReal ratio; 60 PetscInt dim, Nv, v, cStart, cEnd, c; 61 PetscBool adapt = PETSC_TRUE; 62 PetscErrorCode ierr; 63 64 PetscFunctionBeginUser; 65 if (!ctx->adapt) PetscFunctionReturn(0); 66 ierr = DMHasLabel(*dm, "Cell Sets", &hasLabel);CHKERRQ(ierr); 67 if (!hasLabel) {ierr = CreateDomainLabel(*dm);CHKERRQ(ierr);} 68 ierr = DMGetDimension(*dm, &dim);CHKERRQ(ierr); 69 ratio = PetscPowRealInt(0.5, dim); 70 /* Get volume constraints */ 71 ierr = DMGetLabel(*dm, "Cell Sets", &label);CHKERRQ(ierr); 72 ierr = DMLabelGetValueIS(label, &vIS);CHKERRQ(ierr); 73 ierr = ISDuplicate(vIS, &valueIS);CHKERRQ(ierr); 74 ierr = ISDestroy(&vIS);CHKERRQ(ierr); 75 /* Sorting ruins the label */ 76 ierr = ISSort(valueIS);CHKERRQ(ierr); 77 ierr = ISGetLocalSize(valueIS, &Nv);CHKERRQ(ierr); 78 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 79 ierr = PetscMalloc1(Nv, &volConst);CHKERRQ(ierr); 80 for (v = 0; v < Nv; ++v) { 81 char opt[128]; 82 83 volConst[v] = PETSC_MAX_REAL; 84 ierr = PetscSNPrintf(opt, 128, "-volume_constraint_%d", (int) values[v]);CHKERRQ(ierr); 85 ierr = PetscOptionsGetReal(NULL, NULL, opt, &volConst[v], NULL);CHKERRQ(ierr); 86 } 87 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 88 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 89 /* Adapt mesh iteratively */ 90 while (adapt) { 91 DM dmAdapt; 92 DMLabel adaptLabel; 93 PetscInt nAdaptLoc[2], nAdapt[2]; 94 95 adapt = PETSC_FALSE; 96 nAdaptLoc[0] = nAdaptLoc[1] = 0; 97 nAdapt[0] = nAdapt[1] = 0; 98 /* Adaptation is not preserving the domain label */ 99 ierr = DMHasLabel(dmCur, "Cell Sets", &hasLabel);CHKERRQ(ierr); 100 if (!hasLabel) {ierr = CreateDomainLabel(dmCur);CHKERRQ(ierr);} 101 ierr = DMGetLabel(dmCur, "Cell Sets", &label);CHKERRQ(ierr); 102 ierr = DMLabelGetValueIS(label, &vIS);CHKERRQ(ierr); 103 ierr = ISDuplicate(vIS, &valueIS);CHKERRQ(ierr); 104 ierr = ISDestroy(&vIS);CHKERRQ(ierr); 105 /* Sorting directly the label's value IS would corrupt the label so we duplicate the IS first */ 106 ierr = ISSort(valueIS);CHKERRQ(ierr); 107 ierr = ISGetLocalSize(valueIS, &Nv);CHKERRQ(ierr); 108 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 109 /* Construct adaptation label */ 110 ierr = DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel);CHKERRQ(ierr); 111 ierr = DMPlexGetHeightStratum(dmCur, 0, &cStart, &cEnd);CHKERRQ(ierr); 112 for (c = cStart; c < cEnd; ++c) { 113 PetscReal volume, centroid[3]; 114 PetscInt value, vidx; 115 116 ierr = DMPlexComputeCellGeometryFVM(dmCur, c, &volume, centroid, NULL);CHKERRQ(ierr); 117 ierr = DMLabelGetValue(label, c, &value);CHKERRQ(ierr); 118 if (value < 0) continue; 119 ierr = PetscFindInt(value, Nv, values, &vidx);CHKERRQ(ierr); 120 if (vidx < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Value %D for cell %D does not exist in label", value, c); 121 if (volume > volConst[vidx]) {ierr = DMLabelSetValue(adaptLabel, c, DM_ADAPT_REFINE);CHKERRQ(ierr); ++nAdaptLoc[0];} 122 if (volume < volConst[vidx]*ratio) {ierr = DMLabelSetValue(adaptLabel, c, DM_ADAPT_COARSEN);CHKERRQ(ierr); ++nAdaptLoc[1];} 123 } 124 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 125 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 126 ierr = MPI_Allreduce(&nAdaptLoc, &nAdapt, 2, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) dmCur));CHKERRQ(ierr); 127 if (nAdapt[0]) { 128 ierr = PetscInfo2(dmCur, "Adapted mesh, marking %D cells for refinement, and %D cells for coarsening\n", nAdapt[0], nAdapt[1]);CHKERRQ(ierr); 129 ierr = DMAdaptLabel(dmCur, adaptLabel, &dmAdapt);CHKERRQ(ierr); 130 ierr = DMDestroy(&dmCur); 131 ierr = DMViewFromOptions(dmAdapt, NULL, "-adapt_dm_view");CHKERRQ(ierr); 132 dmCur = dmAdapt; 133 adapt = PETSC_TRUE; 134 } 135 ierr = DMLabelDestroy(&adaptLabel);CHKERRQ(ierr); 136 } 137 ierr = PetscFree(volConst);CHKERRQ(ierr); 138 *dm = dmCur; 139 PetscFunctionReturn(0); 140 } 141 142 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 143 { 144 PetscInt dim; 145 size_t len; 146 PetscErrorCode ierr; 147 148 PetscFunctionBeginUser; 149 /* Create top surface */ 150 ierr = PetscStrlen(user->filename, &len);CHKERRQ(ierr); 151 if (len) {ierr = DMPlexCreateFromFile(comm, user->filename, PETSC_TRUE, dm);CHKERRQ(ierr);} 152 else {ierr = DMPlexCreateBoxMesh(comm, 2, PETSC_TRUE, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);CHKERRQ(ierr);} 153 /* Adapt surface */ 154 ierr = AdaptMesh(dm, user);CHKERRQ(ierr); 155 /* Extrude surface to get volume mesh */ 156 ierr = DMGetDimension(*dm, &dim);CHKERRQ(ierr); 157 if (dim < 3 && user->extrude) { 158 DM edm; 159 160 ierr = PetscObjectSetOptionsPrefix((PetscObject) *dm, "srf_");CHKERRQ(ierr); 161 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 162 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 163 ierr = PetscObjectSetOptionsPrefix((PetscObject) *dm, NULL);CHKERRQ(ierr); 164 ierr = DMPlexExtrude(*dm, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_TRUE, NULL, PETSC_TRUE, &edm);CHKERRQ(ierr); 165 ierr = DMDestroy(dm);CHKERRQ(ierr); 166 *dm = edm; 167 } 168 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); 169 ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr); 170 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 171 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 172 PetscFunctionReturn(0); 173 } 174 175 int main(int argc, char **argv) 176 { 177 DM dm; 178 AppCtx user; 179 PetscErrorCode ierr; 180 181 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 182 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 183 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 184 ierr = DMDestroy(&dm);CHKERRQ(ierr); 185 ierr = PetscFinalize(); 186 return ierr; 187 } 188 189 /*TEST 190 191 test: 192 suffix: 0 193 requires: triangle 194 args: -dm_plex_box_dim 2 -dm_plex_box_faces 1,1 -dm_view 195 196 test: # Regularly refine the surface before extrusion 197 suffix: 1 198 requires: triangle 199 args: -dm_plex_box_dim 2 -srf_dm_refine 2 -dm_view 200 201 test: # Parallel run 202 suffix: 2 203 requires: triangle 204 nsize: 5 205 args: -dm_plex_box_dim 2 -srf_dm_refine 3 -petscpartitioner_type simple -dm_distribute -dm_plex_extrude_layers 3 -dm_view 206 207 test: # adaptively refine the surface before extrusion 208 suffix: 3 209 requires: triangle 210 args: -dm_plex_box_dim 2 -dm_plex_box_faces 5,5 -adapt -volume_constraint_1 0.01 -volume_constraint_2 0.000625 -dm_plex_extrude_layers 10 211 212 TEST*/ 213