1 static char help[] = "Test metric utils in the uniform, isotropic case.\n\n"; 2 3 #include <petscdmplex.h> 4 5 static PetscErrorCode bowl(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 6 { 7 PetscInt d; 8 9 *u = 0.0; 10 for (d = 0; d < dim; d++) *u += 0.5*(x[d] - 0.5)*(x[d] - 0.5); 11 12 return 0; 13 } 14 15 static PetscErrorCode CreateIndicator(DM dm, Vec *indicator, DM *dmIndi) 16 { 17 MPI_Comm comm; 18 PetscErrorCode ierr; 19 PetscFE fe; 20 PetscInt dim; 21 22 PetscFunctionBeginUser; 23 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 24 ierr = DMClone(dm, dmIndi);CHKERRQ(ierr); 25 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 26 ierr = PetscFECreateLagrange(comm, dim, 1, PETSC_TRUE, 1, PETSC_DETERMINE, &fe);CHKERRQ(ierr); 27 ierr = DMSetField(*dmIndi, 0, NULL, (PetscObject)fe);CHKERRQ(ierr); 28 ierr = DMCreateDS(*dmIndi);CHKERRQ(ierr); 29 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 30 ierr = DMCreateLocalVector(*dmIndi, indicator);CHKERRQ(ierr); 31 PetscFunctionReturn(0); 32 } 33 34 int main(int argc, char **argv) { 35 DM dm, dmDist, dmAdapt; 36 DMLabel bdLabel = NULL, rgLabel = NULL; 37 MPI_Comm comm; 38 PetscBool uniform = PETSC_FALSE, isotropic = PETSC_FALSE, noTagging = PETSC_FALSE; 39 PetscErrorCode ierr; 40 PetscInt *faces, dim = 3, d; 41 PetscReal scaling = 1.0; 42 Vec metric; 43 44 /* Set up */ 45 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 46 comm = PETSC_COMM_WORLD; 47 ierr = PetscOptionsBegin(comm, "", "Mesh adaptation options", "DMPLEX");CHKERRQ(ierr); 48 ierr = PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex60.c", dim, &dim, NULL, 2, 3);CHKERRQ(ierr); 49 ierr = PetscOptionsBool("-noTagging", "Should tag preservation testing be turned off?", "ex60.c", noTagging, &noTagging, NULL);CHKERRQ(ierr); 50 ierr = PetscOptionsEnd(); 51 52 /* Create box mesh */ 53 ierr = PetscMalloc1(dim, &faces);CHKERRQ(ierr); 54 for (d = 0; d < dim; ++d) faces[d] = 4; 55 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_TRUE, faces, NULL, NULL, NULL, PETSC_TRUE, &dm);CHKERRQ(ierr); 56 ierr = PetscFree(faces);CHKERRQ(ierr); 57 ierr = DMSetFromOptions(dm);CHKERRQ(ierr); 58 59 /* Distribute mesh over processes */ 60 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist);CHKERRQ(ierr); 61 if (dmDist) { 62 ierr = DMDestroy(&dm);CHKERRQ(ierr); 63 dm = dmDist; 64 } 65 ierr = PetscObjectSetName((PetscObject) dm, "DM_init");CHKERRQ(ierr); 66 ierr = DMViewFromOptions(dm, NULL, "-initial_mesh_view");CHKERRQ(ierr); 67 68 /* Set tags to be preserved */ 69 if (!noTagging) { 70 DM cdm; 71 PetscInt cStart, cEnd, c, fStart, fEnd, f, vStart, vEnd; 72 const PetscScalar *coords; 73 Vec coordinates; 74 75 /* Cell tags */ 76 ierr = DMCreateLabel(dm, "Cell Sets");CHKERRQ(ierr); 77 ierr = DMGetLabel(dm, "Cell Sets", &rgLabel);CHKERRQ(ierr); 78 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 79 for (c = cStart; c < cEnd; ++c) { 80 PetscReal centroid[3], volume, x; 81 82 ierr = DMPlexComputeCellGeometryFVM(dm, c, &volume, centroid, NULL);CHKERRQ(ierr); 83 x = centroid[0]; 84 if (x < 0.5) { ierr = DMLabelSetValue(rgLabel, c, 3);CHKERRQ(ierr); } 85 else { ierr = DMLabelSetValue(rgLabel, c, 4);CHKERRQ(ierr); } 86 } 87 88 /* Face tags */ 89 ierr = DMCreateLabel(dm, "Face Sets");CHKERRQ(ierr); 90 ierr = DMGetLabel(dm, "Face Sets", &bdLabel);CHKERRQ(ierr); 91 ierr = DMPlexMarkBoundaryFaces(dm, 1, bdLabel);CHKERRQ(ierr); 92 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 93 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 94 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 95 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 96 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr); 97 for (f = fStart; f < fEnd; ++f) { 98 PetscBool flg = PETSC_TRUE; 99 PetscInt *closure = NULL, closureSize, cl; 100 PetscReal eps = 1.0e-08; 101 102 ierr = DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 103 for (cl = 0; cl < closureSize*2; cl += 2) { 104 PetscInt off = closure[cl]; 105 PetscReal *x; 106 107 if ((off < vStart) || (off >= vEnd)) continue; 108 ierr = DMPlexPointLocalRead(cdm, off, coords, &x);CHKERRQ(ierr); 109 if ((x[0] < 0.5 - eps) || (x[0] > 0.5 + eps)) flg = PETSC_FALSE; 110 } 111 if (flg) { ierr = DMLabelSetValue(bdLabel, f, 2);CHKERRQ(ierr); } 112 ierr = DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 113 } 114 ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr); 115 } 116 117 /* Construct metric */ 118 ierr = DMPlexMetricSetFromOptions(dm);CHKERRQ(ierr); 119 ierr = DMPlexMetricIsUniform(dm, &uniform);CHKERRQ(ierr); 120 ierr = DMPlexMetricIsIsotropic(dm, &isotropic);CHKERRQ(ierr); 121 if (uniform) { 122 ierr = DMPlexMetricCreateUniform(dm, 0, scaling, &metric);CHKERRQ(ierr); 123 } 124 else { 125 DM dmIndi; 126 Vec indicator; 127 128 /* Construct "error indicator" */ 129 ierr = CreateIndicator(dm, &indicator, &dmIndi);CHKERRQ(ierr); 130 if (isotropic) { 131 132 /* Isotropic case: just specify unity */ 133 ierr = VecSet(indicator, scaling);CHKERRQ(ierr); 134 ierr = DMPlexMetricCreateIsotropic(dm, 0, indicator, &metric);CHKERRQ(ierr); 135 136 } else { 137 PetscFE fe; 138 139 /* 'Anisotropic' case: approximate the identity by recovering the Hessian of a parabola */ 140 DM dmGrad; 141 PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar*, void*) = {bowl}; 142 Vec gradient; 143 144 /* Project the parabola into P1 space */ 145 ierr = DMProjectFunctionLocal(dmIndi, 0.0, funcs, NULL, INSERT_ALL_VALUES, indicator);CHKERRQ(ierr); 146 147 /* Approximate the gradient */ 148 ierr = DMClone(dmIndi, &dmGrad);CHKERRQ(ierr); 149 ierr = PetscFECreateLagrange(comm, dim, dim, PETSC_TRUE, 1, PETSC_DETERMINE, &fe);CHKERRQ(ierr); 150 ierr = DMSetField(dmGrad, 0, NULL, (PetscObject)fe);CHKERRQ(ierr); 151 ierr = DMCreateDS(dmGrad);CHKERRQ(ierr); 152 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 153 ierr = DMCreateLocalVector(dmGrad, &gradient);CHKERRQ(ierr); 154 ierr = DMPlexComputeGradientClementInterpolant(dmIndi, indicator, gradient);CHKERRQ(ierr); 155 ierr = VecViewFromOptions(gradient, NULL, "-adapt_gradient_view");CHKERRQ(ierr); 156 157 /* Approximate the Hessian */ 158 ierr = DMPlexMetricCreate(dm, 0, &metric);CHKERRQ(ierr); 159 ierr = DMPlexComputeGradientClementInterpolant(dmGrad, gradient, metric);CHKERRQ(ierr); 160 ierr = VecViewFromOptions(metric, NULL, "-adapt_hessian_view");CHKERRQ(ierr); 161 ierr = VecDestroy(&gradient);CHKERRQ(ierr); 162 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr); 163 } 164 ierr = VecDestroy(&indicator);CHKERRQ(ierr); 165 ierr = DMDestroy(&dmIndi);CHKERRQ(ierr); 166 } 167 168 /* Test metric routines */ 169 { 170 DM dmDet; 171 PetscReal errornorm, norm, tol = 1.0e-10, weights[2] = {0.8, 0.2}; 172 Vec metric1, metric2, metricComb, determinant; 173 Vec metrics[2]; 174 175 ierr = VecDuplicate(metric, &metric1);CHKERRQ(ierr); 176 ierr = VecSet(metric1, 0);CHKERRQ(ierr); 177 ierr = VecAXPY(metric1, 0.625, metric);CHKERRQ(ierr); 178 ierr = VecDuplicate(metric, &metric2);CHKERRQ(ierr); 179 ierr = VecSet(metric2, 0);CHKERRQ(ierr); 180 ierr = VecAXPY(metric2, 2.5, metric);CHKERRQ(ierr); 181 metrics[0] = metric1; 182 metrics[1] = metric2; 183 184 /* Test metric average */ 185 ierr = DMPlexMetricAverage(dm, 2, weights, metrics, &metricComb);CHKERRQ(ierr); 186 ierr = VecAXPY(metricComb, -1, metric);CHKERRQ(ierr); 187 ierr = VecNorm(metric, NORM_2, &norm);CHKERRQ(ierr); 188 ierr = VecNorm(metricComb, NORM_2, &errornorm);CHKERRQ(ierr); 189 errornorm /= norm; 190 ierr = PetscPrintf(comm, "Metric average L2 error: %.4f%%\n", 100*errornorm);CHKERRQ(ierr); 191 if (errornorm > tol) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Metric average test failed"); 192 ierr = VecDestroy(&metricComb);CHKERRQ(ierr); 193 194 /* Test metric intersection */ 195 if (isotropic) { 196 ierr = DMPlexMetricIntersection(dm, 2, metrics, &metricComb);CHKERRQ(ierr); 197 ierr = VecAXPY(metricComb, -1, metric1);CHKERRQ(ierr); 198 ierr = VecNorm(metricComb, NORM_2, &errornorm);CHKERRQ(ierr); 199 errornorm /= norm; 200 ierr = PetscPrintf(comm, "Metric intersection L2 error: %.4f%%\n", 100*errornorm);CHKERRQ(ierr); 201 if (errornorm > tol) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Metric intersection test failed"); 202 } 203 ierr = VecDestroy(&metric1);CHKERRQ(ierr); 204 ierr = VecDestroy(&metric2);CHKERRQ(ierr); 205 ierr = VecDestroy(&metricComb);CHKERRQ(ierr); 206 207 /* Test metric SPD enforcement */ 208 ierr = DMPlexMetricEnforceSPD(dm, metric, PETSC_TRUE, PETSC_TRUE, &metric1, &determinant);CHKERRQ(ierr); 209 if (isotropic) { 210 Vec err; 211 212 ierr = VecDuplicate(determinant, &err);CHKERRQ(ierr); 213 ierr = VecSet(err, 1.0);CHKERRQ(ierr); 214 ierr = VecNorm(err, NORM_2, &norm);CHKERRQ(ierr); 215 ierr = VecAXPY(err, -1, determinant);CHKERRQ(ierr); 216 ierr = VecNorm(err, NORM_2, &errornorm);CHKERRQ(ierr); 217 ierr = VecDestroy(&err);CHKERRQ(ierr); 218 errornorm /= norm; 219 ierr = PetscPrintf(comm, "Metric determinant L2 error: %.4f%%\n", 100*errornorm);CHKERRQ(ierr); 220 if (errornorm > tol) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Determinant is not unit"); 221 ierr = VecAXPY(metric1, -1, metric);CHKERRQ(ierr); 222 ierr = VecNorm(metric1, NORM_2, &errornorm);CHKERRQ(ierr); 223 errornorm /= norm; 224 ierr = PetscPrintf(comm, "Metric SPD enforcement L2 error: %.4f%%\n", 100*errornorm);CHKERRQ(ierr); 225 if (errornorm > tol) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Metric SPD enforcement test failed"); 226 } 227 ierr = VecDestroy(&metric1);CHKERRQ(ierr); 228 ierr = VecGetDM(determinant, &dmDet);CHKERRQ(ierr); 229 ierr = VecDestroy(&determinant);CHKERRQ(ierr); 230 ierr = DMDestroy(&dmDet);CHKERRQ(ierr); 231 232 /* Test metric normalization */ 233 ierr = DMPlexMetricNormalize(dm, metric, PETSC_TRUE, PETSC_TRUE, &metric1);CHKERRQ(ierr); 234 if (isotropic) { 235 PetscReal target; 236 237 ierr = DMPlexMetricGetTargetComplexity(dm, &target);CHKERRQ(ierr); 238 scaling = PetscPowReal(target, 2.0/dim); 239 if (uniform) { 240 ierr = DMPlexMetricCreateUniform(dm, 0, scaling, &metric2);CHKERRQ(ierr); 241 } else { 242 DM dmIndi; 243 Vec indicator; 244 245 ierr = CreateIndicator(dm, &indicator, &dmIndi);CHKERRQ(ierr); 246 ierr = VecSet(indicator, scaling);CHKERRQ(ierr); 247 ierr = DMPlexMetricCreateIsotropic(dm, 0, indicator, &metric2);CHKERRQ(ierr); 248 ierr = DMDestroy(&dmIndi);CHKERRQ(ierr); 249 ierr = VecDestroy(&indicator);CHKERRQ(ierr); 250 } 251 ierr = VecAXPY(metric2, -1, metric1);CHKERRQ(ierr); 252 ierr = VecNorm(metric2, NORM_2, &errornorm);CHKERRQ(ierr); 253 errornorm /= norm; 254 ierr = PetscPrintf(comm, "Metric normalization L2 error: %.4f%%\n", 100*errornorm);CHKERRQ(ierr); 255 if (errornorm > tol) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Metric normalization test failed"); 256 } 257 ierr = VecCopy(metric1, metric);CHKERRQ(ierr); 258 ierr = VecDestroy(&metric2);CHKERRQ(ierr); 259 ierr = VecDestroy(&metric1);CHKERRQ(ierr); 260 } 261 262 /* Adapt the mesh */ 263 ierr = DMAdaptMetric(dm, metric, bdLabel, rgLabel, &dmAdapt);CHKERRQ(ierr); 264 ierr = DMDestroy(&dm);CHKERRQ(ierr); 265 ierr = PetscObjectSetName((PetscObject) dmAdapt, "DM_adapted");CHKERRQ(ierr); 266 ierr = VecDestroy(&metric);CHKERRQ(ierr); 267 ierr = DMViewFromOptions(dmAdapt, NULL, "-adapted_mesh_view");CHKERRQ(ierr); 268 269 /* Test tag preservation */ 270 if (!noTagging) { 271 PetscBool hasTag; 272 PetscInt size; 273 274 ierr = DMGetLabel(dmAdapt, "Face Sets", &bdLabel);CHKERRQ(ierr); 275 ierr = DMLabelHasStratum(bdLabel, 1, &hasTag);CHKERRQ(ierr); 276 if (!hasTag) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Adapted mesh does not have face tag 1"); 277 ierr = DMLabelHasStratum(bdLabel, 2, &hasTag);CHKERRQ(ierr); 278 if (!hasTag) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Adapted mesh does not have face tag 2"); 279 ierr = DMLabelGetNumValues(bdLabel, &size);CHKERRQ(ierr); 280 if (size != 2) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Adapted mesh has the wrong number of face tags (got %d, expected 2)", size); 281 282 ierr = DMGetLabel(dmAdapt, "Cell Sets", &rgLabel);CHKERRQ(ierr); 283 ierr = DMLabelHasStratum(rgLabel, 3, &hasTag);CHKERRQ(ierr); 284 if (!hasTag) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Adapted mesh does not have cell tag 3"); 285 ierr = DMLabelHasStratum(rgLabel, 4, &hasTag);CHKERRQ(ierr); 286 if (!hasTag) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Adapted mesh does not have cell tag 4"); 287 ierr = DMLabelGetNumValues(rgLabel, &size);CHKERRQ(ierr); 288 if (size != 2) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Adapted mesh has the wrong number of cell tags (got %d, expected 2)", size); 289 } 290 291 /* Clean up */ 292 ierr = DMDestroy(&dmAdapt);CHKERRQ(ierr); 293 ierr = PetscFinalize(); 294 return 0; 295 } 296 297 /*TEST 298 299 testset: 300 requires: pragmatic 301 args: -dm_plex_metric_target_complexity 100 -dm_adaptor pragmatic -noTagging -dim 2 302 303 test: 304 suffix: uniform_2d_pragmatic 305 args: -dm_plex_metric_uniform 306 test: 307 suffix: iso_2d_pragmatic 308 args: -dm_plex_metric_isotropic 309 test: 310 suffix: hessian_2d_pragmatic 311 312 testset: 313 requires: pragmatic tetgen 314 args: -dm_plex_metric_target_complexity 100 -dm_adaptor pragmatic -noTagging -dim 3 315 316 test: 317 suffix: uniform_3d_pragmatic 318 args: -dm_plex_metric_uniform -noTagging 319 test: 320 suffix: iso_3d_pragmatic 321 args: -dm_plex_metric_isotropic -noTagging 322 test: 323 suffix: hessian_3d_pragmatic 324 325 testset: 326 requires: mmg 327 args: -dm_plex_metric_target_complexity 100 -dm_adaptor mmg -dim 2 328 329 test: 330 suffix: uniform_2d_mmg 331 args: -dm_plex_metric_uniform 332 test: 333 suffix: iso_2d_mmg 334 args: -dm_plex_metric_isotropic 335 test: 336 suffix: hessian_2d_mmg 337 338 testset: 339 requires: mmg tetgen 340 args: -dm_plex_metric_target_complexity 100 -dm_adaptor mmg -dim 3 341 342 test: 343 suffix: uniform_3d_mmg 344 args: -dm_plex_metric_uniform 345 test: 346 suffix: iso_3d_mmg 347 args: -dm_plex_metric_isotropic 348 test: 349 suffix: hessian_3d_mmg 350 351 testset: 352 requires: parmmg tetgen 353 nsize: 2 354 args: -dm_plex_metric_target_complexity 100 -dm_adaptor parmmg -dim 3 355 356 test: 357 suffix: uniform_3d_parmmg 358 args: -dm_plex_metric_uniform 359 test: 360 suffix: iso_3d_parmmg 361 args: -dm_plex_metric_isotropic 362 test: 363 suffix: hessian_3d_parmmg 364 365 TEST*/ 366