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