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