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