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