Actual source code: ex60.c

  1: static char help[] = "Test metric utils in the uniform, isotropic case.\n\n";

  3: #include <petscdmplex.h>

  5: static PetscErrorCode bowl(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
  6: {
  7:   PetscInt d;

  9:   *u = 0.0;
 10:   for (d = 0; d < dim; d++) *u += 0.5 * (x[d] - 0.5) * (x[d] - 0.5);

 12:   return PETSC_SUCCESS;
 13: }

 15: static PetscErrorCode CreateIndicator(DM dm, Vec *indicator, DM *dmIndi)
 16: {
 17:   MPI_Comm comm;
 18:   PetscFE  fe;
 19:   PetscInt dim;

 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(PETSC_SUCCESS);
 31: }

 33: int main(int argc, char **argv)
 34: {
 35:   DM        dm, dmAdapt;
 36:   DMLabel   bdLabel = NULL, rgLabel = NULL;
 37:   MPI_Comm  comm;
 38:   PetscBool uniform = PETSC_FALSE, isotropic = PETSC_FALSE, noTagging = PETSC_FALSE;
 39:   PetscInt  dim;
 40:   PetscReal scaling = 1.0;
 41:   Vec       metric;

 43:   /* Set up */
 44:   PetscFunctionBeginUser;
 45:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 46:   comm = PETSC_COMM_WORLD;
 47:   PetscOptionsBegin(comm, "", "Mesh adaptation options", "DMPLEX");
 48:   PetscCall(PetscOptionsBool("-noTagging", "Should tag preservation testing be turned off?", "ex60.c", noTagging, &noTagging, NULL));
 49:   PetscOptionsEnd();

 51:   /* Create box mesh */
 52:   PetscCall(DMCreate(comm, &dm));
 53:   PetscCall(DMSetType(dm, DMPLEX));
 54:   PetscCall(DMSetFromOptions(dm));
 55:   PetscCall(PetscObjectSetName((PetscObject)dm, "DM_init"));
 56:   PetscCall(DMViewFromOptions(dm, NULL, "-initial_mesh_view"));
 57:   PetscCall(DMGetDimension(dm, &dim));

 59:   /* Set tags to be preserved */
 60:   if (!noTagging) {
 61:     DM                 cdm;
 62:     PetscInt           cStart, cEnd, c, fStart, fEnd, f, vStart, vEnd;
 63:     const PetscScalar *coords;
 64:     Vec                coordinates;

 66:     /* Cell tags */
 67:     PetscCall(DMGetCoordinatesLocalSetUp(dm));
 68:     PetscCall(DMCreateLabel(dm, "Cell Sets"));
 69:     PetscCall(DMGetLabel(dm, "Cell Sets", &rgLabel));
 70:     PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
 71:     for (c = cStart; c < cEnd; ++c) {
 72:       PetscReal centroid[3], volume, x;

 74:       PetscCall(DMPlexComputeCellGeometryFVM(dm, c, &volume, centroid, NULL));
 75:       x = centroid[0];
 76:       if (x < 0.5) PetscCall(DMLabelSetValue(rgLabel, c, 3));
 77:       else PetscCall(DMLabelSetValue(rgLabel, c, 4));
 78:     }

 80:     /* Face tags */
 81:     PetscCall(DMCreateLabel(dm, "Face Sets"));
 82:     PetscCall(DMGetLabel(dm, "Face Sets", &bdLabel));
 83:     PetscCall(DMPlexMarkBoundaryFaces(dm, 1, bdLabel));
 84:     PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
 85:     PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
 86:     PetscCall(DMGetCoordinateDM(dm, &cdm));
 87:     PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
 88:     PetscCall(VecGetArrayRead(coordinates, &coords));
 89:     for (f = fStart; f < fEnd; ++f) {
 90:       PetscBool flg     = PETSC_TRUE;
 91:       PetscInt *closure = NULL, closureSize, cl;
 92:       PetscReal eps     = 1.0e-08;

 94:       PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure));
 95:       for (cl = 0; cl < closureSize * 2; cl += 2) {
 96:         PetscInt   off = closure[cl];
 97:         PetscReal *x;

 99:         if ((off < vStart) || (off >= vEnd)) continue;
100:         PetscCall(DMPlexPointLocalRead(cdm, off, coords, &x));
101:         if ((x[0] < 0.5 - eps) || (x[0] > 0.5 + eps)) flg = PETSC_FALSE;
102:       }
103:       if (flg) PetscCall(DMLabelSetValue(bdLabel, f, 2));
104:       PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure));
105:     }
106:     PetscCall(VecRestoreArrayRead(coordinates, &coords));
107:   }

109:   /* Construct metric */
110:   PetscCall(DMPlexMetricSetFromOptions(dm));
111:   PetscCall(DMPlexMetricIsUniform(dm, &uniform));
112:   PetscCall(DMPlexMetricIsIsotropic(dm, &isotropic));
113:   if (uniform) {
114:     PetscCall(DMPlexMetricCreateUniform(dm, 0, scaling, &metric));
115:   } else {
116:     DM  dmIndi;
117:     Vec indicator;

119:     /* Construct "error indicator" */
120:     PetscCall(CreateIndicator(dm, &indicator, &dmIndi));
121:     if (isotropic) {
122:       /* Isotropic case: just specify unity */
123:       PetscCall(VecSet(indicator, scaling));
124:       PetscCall(DMPlexMetricCreateIsotropic(dm, 0, indicator, &metric));

126:     } else {
127:       PetscFE fe;

129:       /* 'Anisotropic' case: approximate the identity by recovering the Hessian of a parabola */
130:       DM dmGrad;
131:       PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {bowl};
132:       Vec gradient;

134:       /* Project the parabola into P1 space */
135:       PetscCall(DMProjectFunctionLocal(dmIndi, 0.0, funcs, NULL, INSERT_ALL_VALUES, indicator));

137:       /* Approximate the gradient */
138:       PetscCall(DMClone(dmIndi, &dmGrad));
139:       PetscCall(PetscFECreateLagrange(comm, dim, dim, PETSC_TRUE, 1, PETSC_DETERMINE, &fe));
140:       PetscCall(DMSetField(dmGrad, 0, NULL, (PetscObject)fe));
141:       PetscCall(DMCreateDS(dmGrad));
142:       PetscCall(PetscFEDestroy(&fe));
143:       PetscCall(DMCreateLocalVector(dmGrad, &gradient));
144:       PetscCall(DMPlexComputeGradientClementInterpolant(dmIndi, indicator, gradient));
145:       PetscCall(VecViewFromOptions(gradient, NULL, "-adapt_gradient_view"));

147:       /* Approximate the Hessian */
148:       PetscCall(DMPlexMetricCreate(dm, 0, &metric));
149:       PetscCall(DMPlexComputeGradientClementInterpolant(dmGrad, gradient, metric));
150:       PetscCall(VecViewFromOptions(metric, NULL, "-adapt_hessian_view"));
151:       PetscCall(VecDestroy(&gradient));
152:       PetscCall(DMDestroy(&dmGrad));
153:     }
154:     PetscCall(VecDestroy(&indicator));
155:     PetscCall(DMDestroy(&dmIndi));
156:   }

158:   /* Test metric routines */
159:   {
160:     DM        dmDet;
161:     PetscReal errornorm, norm, tol = 1.0e-10, weights[2] = {0.8, 0.2};
162:     Vec       metric1, metric2, metricComb, determinant;
163:     Vec       metrics[2];

165:     PetscCall(VecDuplicate(metric, &metric1));
166:     PetscCall(VecSet(metric1, 0));
167:     PetscCall(VecAXPY(metric1, 0.625, metric));
168:     PetscCall(VecDuplicate(metric, &metric2));
169:     PetscCall(VecSet(metric2, 0));
170:     PetscCall(VecAXPY(metric2, 2.5, metric));
171:     metrics[0] = metric1;
172:     metrics[1] = metric2;

174:     /* Test metric average */
175:     PetscCall(DMPlexMetricCreate(dm, 0, &metricComb));
176:     PetscCall(DMPlexMetricAverage(dm, 2, weights, metrics, metricComb));
177:     PetscCall(VecAXPY(metricComb, -1, metric));
178:     PetscCall(VecNorm(metric, NORM_2, &norm));
179:     PetscCall(VecNorm(metricComb, NORM_2, &errornorm));
180:     errornorm /= norm;
181:     PetscCall(PetscPrintf(comm, "Metric average L2 error: %.4f%%\n", (double)(100 * errornorm)));
182:     PetscCheck(errornorm < tol, comm, PETSC_ERR_ARG_OUTOFRANGE, "Metric average test failed");

184:     /* Test metric intersection */
185:     PetscCall(DMPlexMetricDeterminantCreate(dm, 0, &determinant, &dmDet));
186:     if (!isotropic) {
187:       PetscCall(DMPlexMetricEnforceSPD(dm, metrics[0], PETSC_FALSE, PETSC_FALSE, metricComb, determinant));
188:       PetscCall(VecCopy(metricComb, metrics[0]));
189:       PetscCall(DMPlexMetricEnforceSPD(dm, metrics[1], PETSC_FALSE, PETSC_FALSE, metricComb, determinant));
190:       PetscCall(VecCopy(metricComb, metrics[1]));
191:     }
192:     PetscCall(DMPlexMetricIntersection(dm, 2, metrics, metricComb));
193:     PetscCall(VecAXPY(metricComb, -1, metric2));
194:     PetscCall(VecNorm(metricComb, NORM_2, &errornorm));
195:     errornorm /= norm;
196:     PetscCall(PetscPrintf(comm, "Metric intersection L2 error: %.4f%%\n", (double)(100 * errornorm)));
197:     PetscCheck(errornorm < tol, comm, PETSC_ERR_ARG_OUTOFRANGE, "Metric intersection test failed");
198:     PetscCall(VecDestroy(&metric2));
199:     PetscCall(VecDestroy(&metricComb));

201:     /* Test metric SPD enforcement */
202:     PetscCall(DMPlexMetricEnforceSPD(dm, metric, PETSC_TRUE, PETSC_TRUE, metric1, determinant));
203:     if (isotropic) {
204:       Vec err;

206:       PetscCall(VecDuplicate(determinant, &err));
207:       PetscCall(VecSet(err, 1.0));
208:       PetscCall(VecNorm(err, NORM_2, &norm));
209:       PetscCall(VecAXPY(err, -1, determinant));
210:       PetscCall(VecNorm(err, NORM_2, &errornorm));
211:       PetscCall(VecDestroy(&err));
212:       errornorm /= norm;
213:       PetscCall(PetscPrintf(comm, "Metric determinant L2 error: %.4f%%\n", (double)(100 * errornorm)));
214:       PetscCheck(errornorm < tol, comm, PETSC_ERR_ARG_OUTOFRANGE, "Determinant is not unit");
215:       PetscCall(VecAXPY(metric1, -1, metric));
216:       PetscCall(VecNorm(metric1, NORM_2, &errornorm));
217:       errornorm /= norm;
218:       PetscCall(PetscPrintf(comm, "Metric SPD enforcement L2 error: %.4f%%\n", (double)(100 * errornorm)));
219:       PetscCheck(errornorm < tol, comm, PETSC_ERR_ARG_OUTOFRANGE, "Metric SPD enforcement test failed");
220:     }

222:     /* Test metric normalization */
223:     PetscCall(DMPlexMetricNormalize(dm, metric, PETSC_TRUE, PETSC_TRUE, metric1, determinant));
224:     if (isotropic) {
225:       PetscReal target;

227:       PetscCall(DMPlexMetricGetTargetComplexity(dm, &target));
228:       scaling = PetscPowReal(target, 2.0 / dim);
229:       if (uniform) {
230:         PetscCall(DMPlexMetricCreateUniform(dm, 0, scaling, &metric2));
231:       } else {
232:         DM  dmIndi;
233:         Vec indicator;

235:         PetscCall(CreateIndicator(dm, &indicator, &dmIndi));
236:         PetscCall(VecSet(indicator, scaling));
237:         PetscCall(DMPlexMetricCreateIsotropic(dm, 0, indicator, &metric2));
238:         PetscCall(DMDestroy(&dmIndi));
239:         PetscCall(VecDestroy(&indicator));
240:       }
241:       PetscCall(VecAXPY(metric2, -1, metric1));
242:       PetscCall(VecNorm(metric2, NORM_2, &errornorm));
243:       errornorm /= norm;
244:       PetscCall(PetscPrintf(comm, "Metric normalization L2 error: %.4f%%\n", (double)(100 * errornorm)));
245:       PetscCheck(errornorm < tol, comm, PETSC_ERR_ARG_OUTOFRANGE, "Metric normalization test failed");
246:     }
247:     PetscCall(VecDestroy(&determinant));
248:     PetscCall(DMDestroy(&dmDet));
249:     PetscCall(VecCopy(metric1, metric));
250:     PetscCall(VecDestroy(&metric2));
251:     PetscCall(VecDestroy(&metric1));
252:   }

254:   /* Adapt the mesh */
255:   PetscCall(DMAdaptMetric(dm, metric, bdLabel, rgLabel, &dmAdapt));
256:   PetscCall(DMDestroy(&dm));
257:   PetscCall(PetscObjectSetName((PetscObject)dmAdapt, "DM_adapted"));
258:   PetscCall(VecDestroy(&metric));
259:   PetscCall(DMViewFromOptions(dmAdapt, NULL, "-adapted_mesh_view"));

261:   /* Test tag preservation */
262:   if (!noTagging) {
263:     PetscBool hasTag;
264:     PetscInt  size;

266:     PetscCall(DMGetLabel(dmAdapt, "Face Sets", &bdLabel));
267:     PetscCall(DMLabelHasStratum(bdLabel, 1, &hasTag));
268:     PetscCheck(hasTag, comm, PETSC_ERR_ARG_OUTOFRANGE, "Adapted mesh does not have face tag 1");
269:     PetscCall(DMLabelHasStratum(bdLabel, 2, &hasTag));
270:     PetscCheck(hasTag, comm, PETSC_ERR_ARG_OUTOFRANGE, "Adapted mesh does not have face tag 2");
271:     PetscCall(DMLabelGetNumValues(bdLabel, &size));
272:     PetscCheck(size == 2, comm, PETSC_ERR_ARG_OUTOFRANGE, "Adapted mesh has the wrong number of face tags (got %" PetscInt_FMT ", expected 2)", size);

274:     PetscCall(DMGetLabel(dmAdapt, "Cell Sets", &rgLabel));
275:     PetscCall(DMLabelHasStratum(rgLabel, 3, &hasTag));
276:     PetscCheck(hasTag, comm, PETSC_ERR_ARG_OUTOFRANGE, "Adapted mesh does not have cell tag 3");
277:     PetscCall(DMLabelHasStratum(rgLabel, 4, &hasTag));
278:     PetscCheck(hasTag, comm, PETSC_ERR_ARG_OUTOFRANGE, "Adapted mesh does not have cell tag 4");
279:     PetscCall(DMLabelGetNumValues(rgLabel, &size));
280:     PetscCheck(size == 2, comm, PETSC_ERR_ARG_OUTOFRANGE, "Adapted mesh has the wrong number of cell tags (got %" PetscInt_FMT ", expected 2)", size);
281:   }

283:   /* Clean up */
284:   PetscCall(DMDestroy(&dmAdapt));
285:   PetscCall(PetscFinalize());
286:   return 0;
287: }

289: /*TEST

291:   testset:
292:     requires: pragmatic
293:     args: -dm_plex_box_faces 4,4 -dm_plex_metric_target_complexity 100 -dm_adaptor pragmatic -noTagging

295:     test:
296:       suffix: uniform_2d_pragmatic
297:       args: -dm_plex_metric_uniform
298:     test:
299:       suffix: iso_2d_pragmatic
300:       args: -dm_plex_metric_isotropic
301:     test:
302:       suffix: hessian_2d_pragmatic

304:   testset:
305:     requires: pragmatic tetgen
306:     args: -dm_plex_dim 3 -dm_plex_box_faces 4,4,4 -dm_plex_metric_target_complexity 100 -dm_adaptor pragmatic -noTagging

308:     test:
309:       suffix: uniform_3d_pragmatic
310:       args: -dm_plex_metric_uniform -noTagging
311:     test:
312:       suffix: iso_3d_pragmatic
313:       args: -dm_plex_metric_isotropic -noTagging
314:     test:
315:       suffix: hessian_3d_pragmatic

317:   testset:
318:     requires: mmg
319:     args: -dm_plex_box_faces 4,4 -dm_plex_metric_target_complexity 100 -dm_adaptor mmg

321:     test:
322:       suffix: uniform_2d_mmg
323:       args: -dm_plex_metric_uniform
324:     test:
325:       suffix: iso_2d_mmg
326:       args: -dm_plex_metric_isotropic
327:     test:
328:       suffix: hessian_2d_mmg

330:   testset:
331:     requires: mmg tetgen
332:     args: -dm_plex_dim 3 -dm_plex_box_faces 4,4,4 -dm_plex_metric_target_complexity 100 -dm_adaptor mmg

334:     test:
335:       suffix: uniform_3d_mmg
336:       args: -dm_plex_metric_uniform
337:     test:
338:       suffix: iso_3d_mmg
339:       args: -dm_plex_metric_isotropic
340:     test:
341:       suffix: hessian_3d_mmg

343:   testset:
344:     requires: parmmg tetgen
345:     nsize: 2
346:     args: -dm_plex_dim 3 -dm_plex_box_faces 4,4,4 -dm_plex_metric_target_complexity 100 -dm_adaptor parmmg

348:     test:
349:       suffix: uniform_3d_parmmg
350:       args: -dm_plex_metric_uniform
351:     test:
352:       suffix: iso_3d_parmmg
353:       args: -dm_plex_metric_isotropic
354:     test:
355:       suffix: hessian_3d_parmmg

357: TEST*/