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*/