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