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