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