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