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 0;
13: }
15: static PetscErrorCode CreateIndicator(DM dm, Vec *indicator, DM *dmIndi)
16: {
17: MPI_Comm comm;
18: PetscFE fe;
19: PetscInt dim;
22: PetscObjectGetComm((PetscObject)dm, &comm);
23: DMClone(dm, dmIndi);
24: DMGetDimension(dm, &dim);
25: PetscFECreateLagrange(comm, dim, 1, PETSC_TRUE, 1, PETSC_DETERMINE, &fe);
26: DMSetField(*dmIndi, 0, NULL, (PetscObject)fe);
27: DMCreateDS(*dmIndi);
28: PetscFEDestroy(&fe);
29: DMCreateLocalVector(*dmIndi, indicator);
30: return 0;
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 */
45: PetscInitialize(&argc, &argv, NULL, help);
46: comm = PETSC_COMM_WORLD;
47: PetscOptionsBegin(comm, "", "Mesh adaptation options", "DMPLEX");
48: PetscOptionsBool("-noTagging", "Should tag preservation testing be turned off?", "ex60.c", noTagging, &noTagging, NULL);
49: PetscOptionsEnd();
51: /* Create box mesh */
52: DMCreate(comm, &dm);
53: DMSetType(dm, DMPLEX);
54: DMSetFromOptions(dm);
55: PetscObjectSetName((PetscObject)dm, "DM_init");
56: DMViewFromOptions(dm, NULL, "-initial_mesh_view");
57: 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: DMGetCoordinatesLocalSetUp(dm);
68: DMCreateLabel(dm, "Cell Sets");
69: DMGetLabel(dm, "Cell Sets", &rgLabel);
70: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
71: for (c = cStart; c < cEnd; ++c) {
72: PetscReal centroid[3], volume, x;
74: DMPlexComputeCellGeometryFVM(dm, c, &volume, centroid, NULL);
75: x = centroid[0];
76: if (x < 0.5) DMLabelSetValue(rgLabel, c, 3);
77: else DMLabelSetValue(rgLabel, c, 4);
78: }
80: /* Face tags */
81: DMCreateLabel(dm, "Face Sets");
82: DMGetLabel(dm, "Face Sets", &bdLabel);
83: DMPlexMarkBoundaryFaces(dm, 1, bdLabel);
84: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
85: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
86: DMGetCoordinateDM(dm, &cdm);
87: DMGetCoordinatesLocal(dm, &coordinates);
88: 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: 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: DMPlexPointLocalRead(cdm, off, coords, &x);
101: if ((x[0] < 0.5 - eps) || (x[0] > 0.5 + eps)) flg = PETSC_FALSE;
102: }
103: if (flg) DMLabelSetValue(bdLabel, f, 2);
104: DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);
105: }
106: VecRestoreArrayRead(coordinates, &coords);
107: }
109: /* Construct metric */
110: DMPlexMetricSetFromOptions(dm);
111: DMPlexMetricIsUniform(dm, &uniform);
112: DMPlexMetricIsIsotropic(dm, &isotropic);
113: if (uniform) {
114: DMPlexMetricCreateUniform(dm, 0, scaling, &metric);
115: } else {
116: DM dmIndi;
117: Vec indicator;
119: /* Construct "error indicator" */
120: CreateIndicator(dm, &indicator, &dmIndi);
121: if (isotropic) {
122: /* Isotropic case: just specify unity */
123: VecSet(indicator, scaling);
124: 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: DMProjectFunctionLocal(dmIndi, 0.0, funcs, NULL, INSERT_ALL_VALUES, indicator);
137: /* Approximate the gradient */
138: DMClone(dmIndi, &dmGrad);
139: PetscFECreateLagrange(comm, dim, dim, PETSC_TRUE, 1, PETSC_DETERMINE, &fe);
140: DMSetField(dmGrad, 0, NULL, (PetscObject)fe);
141: DMCreateDS(dmGrad);
142: PetscFEDestroy(&fe);
143: DMCreateLocalVector(dmGrad, &gradient);
144: DMPlexComputeGradientClementInterpolant(dmIndi, indicator, gradient);
145: VecViewFromOptions(gradient, NULL, "-adapt_gradient_view");
147: /* Approximate the Hessian */
148: DMPlexMetricCreate(dm, 0, &metric);
149: DMPlexComputeGradientClementInterpolant(dmGrad, gradient, metric);
150: VecViewFromOptions(metric, NULL, "-adapt_hessian_view");
151: VecDestroy(&gradient);
152: DMDestroy(&dmGrad);
153: }
154: VecDestroy(&indicator);
155: 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: VecDuplicate(metric, &metric1);
166: VecSet(metric1, 0);
167: VecAXPY(metric1, 0.625, metric);
168: VecDuplicate(metric, &metric2);
169: VecSet(metric2, 0);
170: VecAXPY(metric2, 2.5, metric);
171: metrics[0] = metric1;
172: metrics[1] = metric2;
174: /* Test metric average */
175: DMPlexMetricCreate(dm, 0, &metricComb);
176: DMPlexMetricAverage(dm, 2, weights, metrics, metricComb);
177: VecAXPY(metricComb, -1, metric);
178: VecNorm(metric, NORM_2, &norm);
179: VecNorm(metricComb, NORM_2, &errornorm);
180: errornorm /= norm;
181: PetscPrintf(comm, "Metric average L2 error: %.4f%%\n", (double)(100 * errornorm));
184: /* Test metric intersection */
185: DMPlexMetricDeterminantCreate(dm, 0, &determinant, &dmDet);
186: if (!isotropic) {
187: DMPlexMetricEnforceSPD(dm, metrics[0], PETSC_FALSE, PETSC_FALSE, metricComb, determinant);
188: VecCopy(metricComb, metrics[0]);
189: DMPlexMetricEnforceSPD(dm, metrics[1], PETSC_FALSE, PETSC_FALSE, metricComb, determinant);
190: VecCopy(metricComb, metrics[1]);
191: }
192: DMPlexMetricIntersection(dm, 2, metrics, metricComb);
193: VecAXPY(metricComb, -1, metric2);
194: VecNorm(metricComb, NORM_2, &errornorm);
195: errornorm /= norm;
196: PetscPrintf(comm, "Metric intersection L2 error: %.4f%%\n", (double)(100 * errornorm));
198: VecDestroy(&metric2);
199: VecDestroy(&metricComb);
201: /* Test metric SPD enforcement */
202: DMPlexMetricEnforceSPD(dm, metric, PETSC_TRUE, PETSC_TRUE, metric1, determinant);
203: if (isotropic) {
204: Vec err;
206: VecDuplicate(determinant, &err);
207: VecSet(err, 1.0);
208: VecNorm(err, NORM_2, &norm);
209: VecAXPY(err, -1, determinant);
210: VecNorm(err, NORM_2, &errornorm);
211: VecDestroy(&err);
212: errornorm /= norm;
213: PetscPrintf(comm, "Metric determinant L2 error: %.4f%%\n", (double)(100 * errornorm));
215: VecAXPY(metric1, -1, metric);
216: VecNorm(metric1, NORM_2, &errornorm);
217: errornorm /= norm;
218: PetscPrintf(comm, "Metric SPD enforcement L2 error: %.4f%%\n", (double)(100 * errornorm));
220: }
222: /* Test metric normalization */
223: DMPlexMetricNormalize(dm, metric, PETSC_TRUE, PETSC_TRUE, metric1, determinant);
224: if (isotropic) {
225: PetscReal target;
227: DMPlexMetricGetTargetComplexity(dm, &target);
228: scaling = PetscPowReal(target, 2.0 / dim);
229: if (uniform) {
230: DMPlexMetricCreateUniform(dm, 0, scaling, &metric2);
231: } else {
232: DM dmIndi;
233: Vec indicator;
235: CreateIndicator(dm, &indicator, &dmIndi);
236: VecSet(indicator, scaling);
237: DMPlexMetricCreateIsotropic(dm, 0, indicator, &metric2);
238: DMDestroy(&dmIndi);
239: VecDestroy(&indicator);
240: }
241: VecAXPY(metric2, -1, metric1);
242: VecNorm(metric2, NORM_2, &errornorm);
243: errornorm /= norm;
244: PetscPrintf(comm, "Metric normalization L2 error: %.4f%%\n", (double)(100 * errornorm));
246: }
247: VecDestroy(&determinant);
248: DMDestroy(&dmDet);
249: VecCopy(metric1, metric);
250: VecDestroy(&metric2);
251: VecDestroy(&metric1);
252: }
254: /* Adapt the mesh */
255: DMAdaptMetric(dm, metric, bdLabel, rgLabel, &dmAdapt);
256: DMDestroy(&dm);
257: PetscObjectSetName((PetscObject)dmAdapt, "DM_adapted");
258: VecDestroy(&metric);
259: DMViewFromOptions(dmAdapt, NULL, "-adapted_mesh_view");
261: /* Test tag preservation */
262: if (!noTagging) {
263: PetscBool hasTag;
264: PetscInt size;
266: DMGetLabel(dmAdapt, "Face Sets", &bdLabel);
267: DMLabelHasStratum(bdLabel, 1, &hasTag);
269: DMLabelHasStratum(bdLabel, 2, &hasTag);
271: DMLabelGetNumValues(bdLabel, &size);
274: DMGetLabel(dmAdapt, "Cell Sets", &rgLabel);
275: DMLabelHasStratum(rgLabel, 3, &hasTag);
277: DMLabelHasStratum(rgLabel, 4, &hasTag);
279: DMLabelGetNumValues(rgLabel, &size);
281: }
283: /* Clean up */
284: DMDestroy(&dmAdapt);
285: 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*/