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