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