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