1 static char help[] = "Tests mesh adaptation with DMPlex and pragmatic.\n"; 2 3 #include <petsc/private/dmpleximpl.h> 4 5 #include <petscsnes.h> 6 7 typedef struct { 8 PetscInt Nr; /* The number of refinement passes */ 9 PetscInt metOpt; /* Different choices of metric */ 10 PetscReal hmax, hmin; /* Max and min sizes prescribed by the metric */ 11 PetscBool doL2; /* Test L2 projection */ 12 } AppCtx; 13 14 /* 15 Classic hyperbolic sensor function for testing multi-scale anisotropic mesh adaptation: 16 17 f:[-1, 1]x[-1, 1] \to R, 18 f(x, y) = sin(50xy)/100 if |xy| > 2\pi/50 else sin(50xy) 19 20 (mapped to have domain [0,1] x [0,1] in this case). 21 */ 22 static PetscErrorCode sensor(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) { 23 const PetscReal xref = 2. * x[0] - 1.; 24 const PetscReal yref = 2. * x[1] - 1.; 25 const PetscReal xy = xref * yref; 26 27 PetscFunctionBeginUser; 28 u[0] = PetscSinReal(50. * xy); 29 if (PetscAbsReal(xy) > 2. * PETSC_PI / 50.) u[0] *= 0.01; 30 PetscFunctionReturn(0); 31 } 32 33 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) { 34 PetscFunctionBegin; 35 options->Nr = 1; 36 options->metOpt = 1; 37 options->hmin = 0.05; 38 options->hmax = 0.5; 39 options->doL2 = PETSC_FALSE; 40 41 PetscOptionsBegin(comm, "", "Meshing Adaptation Options", "DMPLEX"); 42 PetscCall(PetscOptionsBoundedInt("-Nr", "Numberof refinement passes", "ex19.c", options->Nr, &options->Nr, NULL, 1)); 43 PetscCall(PetscOptionsBoundedInt("-met", "Different choices of metric", "ex19.c", options->metOpt, &options->metOpt, NULL, 0)); 44 PetscCall(PetscOptionsReal("-hmax", "Max size prescribed by the metric", "ex19.c", options->hmax, &options->hmax, NULL)); 45 PetscCall(PetscOptionsReal("-hmin", "Min size prescribed by the metric", "ex19.c", options->hmin, &options->hmin, NULL)); 46 PetscCall(PetscOptionsBool("-do_L2", "Test L2 projection", "ex19.c", options->doL2, &options->doL2, NULL)); 47 PetscOptionsEnd(); 48 49 PetscFunctionReturn(0); 50 } 51 52 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm) { 53 PetscFunctionBegin; 54 PetscCall(DMCreate(comm, dm)); 55 PetscCall(DMSetType(*dm, DMPLEX)); 56 PetscCall(DMSetFromOptions(*dm)); 57 PetscCall(PetscObjectSetName((PetscObject)*dm, "DMinit")); 58 PetscCall(DMViewFromOptions(*dm, NULL, "-init_dm_view")); 59 PetscFunctionReturn(0); 60 } 61 62 static PetscErrorCode ComputeMetricSensor(DM dm, AppCtx *user, Vec *metric) { 63 PetscSimplePointFunc funcs[1] = {sensor}; 64 DM dmSensor, dmGrad, dmHess, dmDet; 65 PetscFE fe; 66 Vec f, g, H, determinant; 67 PetscBool simplex; 68 PetscInt dim; 69 70 PetscFunctionBegin; 71 PetscCall(DMGetDimension(dm, &dim)); 72 PetscCall(DMPlexIsSimplex(dm, &simplex)); 73 74 PetscCall(DMClone(dm, &dmSensor)); 75 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, 1, simplex, 1, -1, &fe)); 76 PetscCall(DMSetField(dmSensor, 0, NULL, (PetscObject)fe)); 77 PetscCall(PetscFEDestroy(&fe)); 78 PetscCall(DMCreateDS(dmSensor)); 79 PetscCall(DMCreateLocalVector(dmSensor, &f)); 80 PetscCall(DMProjectFunctionLocal(dmSensor, 0., funcs, NULL, INSERT_VALUES, f)); 81 PetscCall(VecViewFromOptions(f, NULL, "-sensor_view")); 82 83 // Recover the gradient of the sensor function 84 PetscCall(DMClone(dm, &dmGrad)); 85 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, dim, simplex, 1, -1, &fe)); 86 PetscCall(DMSetField(dmGrad, 0, NULL, (PetscObject)fe)); 87 PetscCall(PetscFEDestroy(&fe)); 88 PetscCall(DMCreateDS(dmGrad)); 89 PetscCall(DMCreateLocalVector(dmGrad, &g)); 90 PetscCall(DMPlexComputeGradientClementInterpolant(dmSensor, f, g)); 91 PetscCall(VecDestroy(&f)); 92 PetscCall(VecViewFromOptions(g, NULL, "-gradient_view")); 93 94 // Recover the Hessian of the sensor function 95 PetscCall(DMClone(dm, &dmHess)); 96 PetscCall(DMPlexMetricCreate(dmHess, 0, &H)); 97 PetscCall(DMPlexComputeGradientClementInterpolant(dmGrad, g, H)); 98 PetscCall(VecDestroy(&g)); 99 PetscCall(VecViewFromOptions(H, NULL, "-hessian_view")); 100 101 // Obtain a metric by Lp normalization 102 PetscCall(DMPlexMetricCreate(dm, 0, metric)); 103 PetscCall(DMPlexMetricDeterminantCreate(dm, 0, &determinant, &dmDet)); 104 PetscCall(DMPlexMetricNormalize(dmHess, H, PETSC_TRUE, PETSC_TRUE, *metric, determinant)); 105 PetscCall(VecDestroy(&determinant)); 106 PetscCall(DMDestroy(&dmDet)); 107 PetscCall(VecDestroy(&H)); 108 PetscCall(DMDestroy(&dmHess)); 109 PetscCall(DMDestroy(&dmGrad)); 110 PetscCall(DMDestroy(&dmSensor)); 111 PetscFunctionReturn(0); 112 } 113 114 static PetscErrorCode ComputeMetric(DM dm, AppCtx *user, Vec *metric) { 115 PetscReal lambda = 1 / (user->hmax * user->hmax); 116 117 PetscFunctionBeginUser; 118 if (user->metOpt == 0) { 119 /* Specify a uniform, isotropic metric */ 120 PetscCall(DMPlexMetricCreateUniform(dm, 0, lambda, metric)); 121 } else if (user->metOpt == 3) { 122 PetscCall(ComputeMetricSensor(dm, user, metric)); 123 } else { 124 DM cdm; 125 Vec coordinates; 126 const PetscScalar *coords; 127 PetscScalar *met; 128 PetscReal h; 129 PetscInt dim, i, j, vStart, vEnd, v; 130 131 PetscCall(DMPlexMetricCreate(dm, 0, metric)); 132 PetscCall(DMGetDimension(dm, &dim)); 133 PetscCall(DMGetCoordinateDM(dm, &cdm)); 134 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 135 PetscCall(VecGetArrayRead(coordinates, &coords)); 136 PetscCall(VecGetArray(*metric, &met)); 137 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 138 for (v = vStart; v < vEnd; ++v) { 139 PetscScalar *vcoords; 140 PetscScalar *pmet; 141 142 PetscCall(DMPlexPointLocalRead(cdm, v, coords, &vcoords)); 143 switch (user->metOpt) { 144 case 1: h = user->hmax - (user->hmax - user->hmin) * PetscRealPart(vcoords[0]); break; 145 case 2: h = user->hmax * PetscAbsReal(((PetscReal)1.0) - PetscExpReal(-PetscAbsScalar(vcoords[0] - (PetscReal)0.5))) + user->hmin; break; 146 default: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "metOpt = 0, 1, 2 or 3, cannot be %d", user->metOpt); 147 } 148 PetscCall(DMPlexPointLocalRef(dm, v, met, &pmet)); 149 for (i = 0; i < dim; ++i) { 150 for (j = 0; j < dim; ++j) { 151 if (i == j) { 152 if (i == 0) pmet[i * dim + j] = 1 / (h * h); 153 else pmet[i * dim + j] = lambda; 154 } else pmet[i * dim + j] = 0.0; 155 } 156 } 157 } 158 PetscCall(VecRestoreArray(*metric, &met)); 159 PetscCall(VecRestoreArrayRead(coordinates, &coords)); 160 } 161 PetscFunctionReturn(0); 162 } 163 164 static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) { 165 u[0] = x[0] + x[1]; 166 return 0; 167 } 168 169 static PetscErrorCode TestL2Projection(DM dm, DM dma, AppCtx *user) { 170 PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {linear}; 171 DM dmProj, dmaProj; 172 PetscFE fe; 173 KSP ksp; 174 Mat Interp, mass, mass2; 175 Vec u, ua, scaling, rhs, uproj; 176 PetscReal error; 177 PetscBool simplex; 178 PetscInt dim; 179 180 PetscFunctionBeginUser; 181 PetscCall(DMGetDimension(dm, &dim)); 182 PetscCall(DMPlexIsSimplex(dm, &simplex)); 183 184 PetscCall(DMClone(dm, &dmProj)); 185 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, -1, &fe)); 186 PetscCall(DMSetField(dmProj, 0, NULL, (PetscObject)fe)); 187 PetscCall(PetscFEDestroy(&fe)); 188 PetscCall(DMCreateDS(dmProj)); 189 190 PetscCall(DMClone(dma, &dmaProj)); 191 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, -1, &fe)); 192 PetscCall(DMSetField(dmaProj, 0, NULL, (PetscObject)fe)); 193 PetscCall(PetscFEDestroy(&fe)); 194 PetscCall(DMCreateDS(dmaProj)); 195 196 PetscCall(DMGetGlobalVector(dmProj, &u)); 197 PetscCall(DMGetGlobalVector(dmaProj, &ua)); 198 PetscCall(DMGetGlobalVector(dmaProj, &rhs)); 199 PetscCall(DMGetGlobalVector(dmaProj, &uproj)); 200 201 // Interpolate onto original mesh using dual basis 202 PetscCall(DMProjectFunction(dmProj, 0.0, funcs, NULL, INSERT_VALUES, u)); 203 PetscCall(PetscObjectSetName((PetscObject)u, "Original")); 204 PetscCall(VecViewFromOptions(u, NULL, "-orig_vec_view")); 205 PetscCall(DMComputeL2Diff(dmProj, 0.0, funcs, NULL, u, &error)); 206 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Original L2 Error: %g\n", (double)error)); 207 // Interpolate onto NEW mesh using dual basis 208 PetscCall(DMProjectFunction(dmaProj, 0.0, funcs, NULL, INSERT_VALUES, ua)); 209 PetscCall(PetscObjectSetName((PetscObject)ua, "Adapted")); 210 PetscCall(VecViewFromOptions(ua, NULL, "-adapt_vec_view")); 211 PetscCall(DMComputeL2Diff(dmaProj, 0.0, funcs, NULL, ua, &error)); 212 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Adapted L2 Error: %g\n", (double)error)); 213 // Interpolate between meshes using interpolation matrix 214 PetscCall(DMCreateInterpolation(dmProj, dmaProj, &Interp, &scaling)); 215 PetscCall(MatInterpolate(Interp, u, ua)); 216 PetscCall(MatDestroy(&Interp)); 217 PetscCall(VecDestroy(&scaling)); 218 PetscCall(PetscObjectSetName((PetscObject)ua, "Interpolation")); 219 PetscCall(VecViewFromOptions(ua, NULL, "-interp_vec_view")); 220 PetscCall(DMComputeL2Diff(dmaProj, 0.0, funcs, NULL, ua, &error)); 221 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Interpolated L2 Error: %g\n", (double)error)); 222 // L2 projection 223 PetscCall(DMCreateMassMatrix(dmaProj, dmaProj, &mass)); 224 PetscCall(MatViewFromOptions(mass, NULL, "-mass_mat_view")); 225 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp)); 226 PetscCall(KSPSetOperators(ksp, mass, mass)); 227 PetscCall(KSPSetFromOptions(ksp)); 228 // Compute rhs as M f, could also direclty project the analytic function but we might not have it 229 PetscCall(DMCreateMassMatrix(dmProj, dmaProj, &mass2)); 230 PetscCall(MatMult(mass2, u, rhs)); 231 PetscCall(MatDestroy(&mass2)); 232 PetscCall(KSPSolve(ksp, rhs, uproj)); 233 PetscCall(PetscObjectSetName((PetscObject)uproj, "L_2 Projection")); 234 PetscCall(VecViewFromOptions(uproj, NULL, "-proj_vec_view")); 235 PetscCall(DMComputeL2Diff(dmaProj, 0.0, funcs, NULL, uproj, &error)); 236 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Projected L2 Error: %g\n", (double)error)); 237 PetscCall(KSPDestroy(&ksp)); 238 PetscCall(MatDestroy(&mass)); 239 PetscCall(DMRestoreGlobalVector(dmProj, &u)); 240 PetscCall(DMRestoreGlobalVector(dmaProj, &ua)); 241 PetscCall(DMRestoreGlobalVector(dmaProj, &rhs)); 242 PetscCall(DMRestoreGlobalVector(dmaProj, &uproj)); 243 PetscCall(DMDestroy(&dmProj)); 244 PetscCall(DMDestroy(&dmaProj)); 245 PetscFunctionReturn(0); 246 } 247 248 int main(int argc, char *argv[]) { 249 DM dm; 250 AppCtx user; /* user-defined work context */ 251 MPI_Comm comm; 252 DM dma, odm; 253 Vec metric; 254 PetscInt r; 255 256 PetscFunctionBeginUser; 257 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 258 comm = PETSC_COMM_WORLD; 259 PetscCall(ProcessOptions(comm, &user)); 260 PetscCall(CreateMesh(comm, &dm)); 261 262 odm = dm; 263 PetscCall(DMPlexDistributeOverlap(odm, 1, NULL, &dm)); 264 if (!dm) { 265 dm = odm; 266 } else PetscCall(DMDestroy(&odm)); 267 268 for (r = 0; r < user.Nr; ++r) { 269 DMLabel label; 270 271 PetscCall(ComputeMetric(dm, &user, &metric)); 272 PetscCall(DMGetLabel(dm, "marker", &label)); 273 PetscCall(DMAdaptMetric(dm, metric, label, NULL, &dma)); 274 PetscCall(VecDestroy(&metric)); 275 PetscCall(PetscObjectSetName((PetscObject)dma, "DMadapt")); 276 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dma, "adapt_")); 277 PetscCall(DMViewFromOptions(dma, NULL, "-dm_view")); 278 if (user.doL2) PetscCall(TestL2Projection(dm, dma, &user)); 279 PetscCall(DMDestroy(&dm)); 280 dm = dma; 281 } 282 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dm, "final_")); 283 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 284 PetscCall(DMDestroy(&dm)); 285 PetscCall(PetscFinalize()); 286 return 0; 287 } 288 289 /*TEST 290 291 build: 292 requires: pragmatic 293 294 testset: 295 args: -dm_plex_box_faces 4,4,4 -dm_adaptor pragmatic -met 2 -init_dm_view -adapt_dm_view -dm_adaptor pragmatic 296 297 test: 298 suffix: 2d 299 args: -dm_plex_separate_marker 0 300 test: 301 suffix: 2d_sep 302 args: -dm_plex_separate_marker 1 303 test: 304 suffix: 3d 305 args: -dm_plex_dim 3 306 307 # Pragmatic hangs for simple partitioner 308 testset: 309 requires: parmetis 310 args: -dm_plex_box_faces 2,2 -petscpartitioner_type parmetis -met 2 -init_dm_view -adapt_dm_view -dm_adaptor pragmatic 311 312 test: 313 suffix: 2d_parmetis_np2 314 nsize: 2 315 test: 316 suffix: 2d_parmetis_np4 317 nsize: 4 318 319 test: 320 requires: parmetis 321 suffix: 3d_parmetis_met0 322 nsize: 2 323 args: -dm_plex_dim 3 -dm_plex_box_faces 9,9,9 -dm_adaptor pragmatic -petscpartitioner_type parmetis \ 324 -met 0 -hmin 0.01 -hmax 0.03 -init_dm_view -adapt_dm_view -dm_adaptor pragmatic 325 test: 326 requires: parmetis 327 suffix: 3d_parmetis_met2 328 nsize: 2 329 args: -dm_plex_box_faces 19,19 -dm_adaptor pragmatic -petscpartitioner_type parmetis \ 330 -met 2 -hmax 0.5 -hmin 0.001 -init_dm_view -adapt_dm_view -dm_adaptor pragmatic 331 test: 332 suffix: proj2 333 args: -dm_plex_box_faces 2,2 -dm_plex_hash_location -dm_adaptor pragmatic -init_dm_view -adapt_dm_view -do_L2 \ 334 -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu -dm_adaptor pragmatic 335 test: 336 suffix: proj4 337 args: -dm_plex_box_faces 4,4 -dm_plex_hash_location -dm_adaptor pragmatic -init_dm_view -adapt_dm_view -do_L2 \ 338 -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu -dm_adaptor pragmatic 339 340 test: 341 suffix: 2d_met3 342 args: -dm_plex_box_faces 9,9 -met 3 -dm_adaptor pragmatic -init_dm_view -adapt_dm_view \ 343 -dm_plex_metric_h_min 1.e-10 -dm_plex_metric_h_max 1.0e-01 -dm_plex_metric_a_max 1.0e+05 -dm_plex_metric_p 1.0 \ 344 -dm_plex_metric_target_complexity 10000.0 -dm_adaptor pragmatic 345 346 TEST*/ 347