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