1 static char help[] = "Test metric utils in the uniform, isotropic case.\n\n";
2
3 #include <petscdmplex.h>
4
bowl(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)5 static PetscErrorCode bowl(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx 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 PETSC_SUCCESS;
13 }
14
CreateIndicator(DM dm,Vec * indicator,DM * dmIndi)15 static PetscErrorCode CreateIndicator(DM dm, Vec *indicator, DM *dmIndi)
16 {
17 MPI_Comm comm;
18 PetscFE fe;
19 PetscInt dim;
20
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 }
32
main(int argc,char ** argv)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;
42
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();
50
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));
58
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;
65
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;
73
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 }
79
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;
93
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;
98
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 }
108
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;
118
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));
125
126 } else {
127 PetscFE fe;
128
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;
133
134 /* Project the parabola into P1 space */
135 PetscCall(DMProjectFunctionLocal(dmIndi, 0.0, funcs, NULL, INSERT_ALL_VALUES, indicator));
136
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"));
146
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 }
157
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];
164
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;
173
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");
183
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));
200
201 /* Test metric SPD enforcement */
202 PetscCall(DMPlexMetricEnforceSPD(dm, metric, PETSC_TRUE, PETSC_TRUE, metric1, determinant));
203 if (isotropic) {
204 Vec err;
205
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 }
221
222 /* Test metric normalization */
223 PetscCall(DMPlexMetricNormalize(dm, metric, PETSC_TRUE, PETSC_TRUE, metric1, determinant));
224 if (isotropic) {
225 PetscReal target;
226
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;
234
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 }
253
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"));
260
261 /* Test tag preservation */
262 if (!noTagging) {
263 PetscBool hasTag;
264 PetscInt size;
265
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);
273
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 }
282
283 /* Clean up */
284 PetscCall(DMDestroy(&dmAdapt));
285 PetscCall(PetscFinalize());
286 return 0;
287 }
288
289 /*TEST
290
291 testset:
292 requires: pragmatic
293 args: -dm_plex_box_faces 4,4 -dm_plex_metric_target_complexity 100 -dm_adaptor pragmatic -noTagging
294
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
303
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
307
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
316
317 testset:
318 requires: mmg
319 args: -dm_plex_box_faces 4,4 -dm_plex_metric_target_complexity 100 -dm_adaptor mmg
320
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
329
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
333
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
342
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
347
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
356
357 TEST*/
358