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