Actual source code: ex19.c
1: static char help[] = "Tests mesh adaptation with DMPlex and pragmatic.\n";
3: #include <petsc/private/dmpleximpl.h>
5: #include <petscsnes.h>
7: typedef struct {
8: PetscInt Nr; /* The number of refinement passes */
9: PetscInt metOpt; /* Different choices of metric */
10: PetscReal hmax, hmin; /* Max and min sizes prescribed by the metric */
11: PetscBool doL2; /* Test L2 projection */
12: } AppCtx;
14: /*
15: Classic hyperbolic sensor function for testing multi-scale anisotropic mesh adaptation:
17: f:[-1, 1]x[-1, 1] \to R,
18: f(x, y) = sin(50xy)/100 if |xy| > 2\pi/50 else sin(50xy)
20: (mapped to have domain [0,1] x [0,1] in this case).
21: */
22: static PetscErrorCode sensor(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
23: {
24: const PetscReal xref = 2. * x[0] - 1.;
25: const PetscReal yref = 2. * x[1] - 1.;
26: const PetscReal xy = xref * yref;
29: u[0] = PetscSinReal(50. * xy);
30: if (PetscAbsReal(xy) > 2. * PETSC_PI / 50.) u[0] *= 0.01;
31: return 0;
32: }
34: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
35: {
36: options->Nr = 1;
37: options->metOpt = 1;
38: options->hmin = 0.05;
39: options->hmax = 0.5;
40: options->doL2 = PETSC_FALSE;
42: PetscOptionsBegin(comm, "", "Meshing Adaptation Options", "DMPLEX");
43: PetscOptionsBoundedInt("-Nr", "Numberof refinement passes", "ex19.c", options->Nr, &options->Nr, NULL, 1);
44: PetscOptionsBoundedInt("-met", "Different choices of metric", "ex19.c", options->metOpt, &options->metOpt, NULL, 0);
45: PetscOptionsReal("-hmax", "Max size prescribed by the metric", "ex19.c", options->hmax, &options->hmax, NULL);
46: PetscOptionsReal("-hmin", "Min size prescribed by the metric", "ex19.c", options->hmin, &options->hmin, NULL);
47: PetscOptionsBool("-do_L2", "Test L2 projection", "ex19.c", options->doL2, &options->doL2, NULL);
48: PetscOptionsEnd();
50: return 0;
51: }
53: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm)
54: {
55: DMCreate(comm, dm);
56: DMSetType(*dm, DMPLEX);
57: DMSetFromOptions(*dm);
58: PetscObjectSetName((PetscObject)*dm, "DMinit");
59: DMViewFromOptions(*dm, NULL, "-init_dm_view");
60: return 0;
61: }
63: static PetscErrorCode ComputeMetricSensor(DM dm, AppCtx *user, Vec *metric)
64: {
65: PetscSimplePointFunc funcs[1] = {sensor};
66: DM dmSensor, dmGrad, dmHess, dmDet;
67: PetscFE fe;
68: Vec f, g, H, determinant;
69: PetscBool simplex;
70: PetscInt dim;
72: DMGetDimension(dm, &dim);
73: DMPlexIsSimplex(dm, &simplex);
75: DMClone(dm, &dmSensor);
76: PetscFECreateLagrange(PETSC_COMM_SELF, dim, 1, simplex, 1, -1, &fe);
77: DMSetField(dmSensor, 0, NULL, (PetscObject)fe);
78: PetscFEDestroy(&fe);
79: DMCreateDS(dmSensor);
80: DMCreateLocalVector(dmSensor, &f);
81: DMProjectFunctionLocal(dmSensor, 0., funcs, NULL, INSERT_VALUES, f);
82: VecViewFromOptions(f, NULL, "-sensor_view");
84: // Recover the gradient of the sensor function
85: DMClone(dm, &dmGrad);
86: PetscFECreateLagrange(PETSC_COMM_SELF, dim, dim, simplex, 1, -1, &fe);
87: DMSetField(dmGrad, 0, NULL, (PetscObject)fe);
88: PetscFEDestroy(&fe);
89: DMCreateDS(dmGrad);
90: DMCreateLocalVector(dmGrad, &g);
91: DMPlexComputeGradientClementInterpolant(dmSensor, f, g);
92: VecDestroy(&f);
93: VecViewFromOptions(g, NULL, "-gradient_view");
95: // Recover the Hessian of the sensor function
96: DMClone(dm, &dmHess);
97: DMPlexMetricCreate(dmHess, 0, &H);
98: DMPlexComputeGradientClementInterpolant(dmGrad, g, H);
99: VecDestroy(&g);
100: VecViewFromOptions(H, NULL, "-hessian_view");
102: // Obtain a metric by Lp normalization
103: DMPlexMetricCreate(dm, 0, metric);
104: DMPlexMetricDeterminantCreate(dm, 0, &determinant, &dmDet);
105: DMPlexMetricNormalize(dmHess, H, PETSC_TRUE, PETSC_TRUE, *metric, determinant);
106: VecDestroy(&determinant);
107: DMDestroy(&dmDet);
108: VecDestroy(&H);
109: DMDestroy(&dmHess);
110: DMDestroy(&dmGrad);
111: DMDestroy(&dmSensor);
112: return 0;
113: }
115: static PetscErrorCode ComputeMetric(DM dm, AppCtx *user, Vec *metric)
116: {
117: PetscReal lambda = 1 / (user->hmax * user->hmax);
120: if (user->metOpt == 0) {
121: /* Specify a uniform, isotropic metric */
122: DMPlexMetricCreateUniform(dm, 0, lambda, metric);
123: } else if (user->metOpt == 3) {
124: ComputeMetricSensor(dm, user, metric);
125: } else {
126: DM cdm;
127: Vec coordinates;
128: const PetscScalar *coords;
129: PetscScalar *met;
130: PetscReal h;
131: PetscInt dim, i, j, vStart, vEnd, v;
133: DMPlexMetricCreate(dm, 0, metric);
134: DMGetDimension(dm, &dim);
135: DMGetCoordinateDM(dm, &cdm);
136: DMGetCoordinatesLocal(dm, &coordinates);
137: VecGetArrayRead(coordinates, &coords);
138: VecGetArray(*metric, &met);
139: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
140: for (v = vStart; v < vEnd; ++v) {
141: PetscScalar *vcoords;
142: PetscScalar *pmet;
144: DMPlexPointLocalRead(cdm, v, coords, &vcoords);
145: switch (user->metOpt) {
146: case 1:
147: h = user->hmax - (user->hmax - user->hmin) * PetscRealPart(vcoords[0]);
148: break;
149: case 2:
150: h = user->hmax * PetscAbsReal(((PetscReal)1.0) - PetscExpReal(-PetscAbsScalar(vcoords[0] - (PetscReal)0.5))) + user->hmin;
151: break;
152: default:
153: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "metOpt = 0, 1, 2 or 3, cannot be %d", user->metOpt);
154: }
155: DMPlexPointLocalRef(dm, v, met, &pmet);
156: for (i = 0; i < dim; ++i) {
157: for (j = 0; j < dim; ++j) {
158: if (i == j) {
159: if (i == 0) pmet[i * dim + j] = 1 / (h * h);
160: else pmet[i * dim + j] = lambda;
161: } else pmet[i * dim + j] = 0.0;
162: }
163: }
164: }
165: VecRestoreArray(*metric, &met);
166: VecRestoreArrayRead(coordinates, &coords);
167: }
168: return 0;
169: }
171: static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
172: {
173: u[0] = x[0] + x[1];
174: return 0;
175: }
177: static PetscErrorCode TestL2Projection(DM dm, DM dma, AppCtx *user)
178: {
179: PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {linear};
180: DM dmProj, dmaProj;
181: PetscFE fe;
182: KSP ksp;
183: Mat Interp, mass, mass2;
184: Vec u, ua, scaling, rhs, uproj;
185: PetscReal error;
186: PetscBool simplex;
187: PetscInt dim;
190: DMGetDimension(dm, &dim);
191: DMPlexIsSimplex(dm, &simplex);
193: DMClone(dm, &dmProj);
194: PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, -1, &fe);
195: DMSetField(dmProj, 0, NULL, (PetscObject)fe);
196: PetscFEDestroy(&fe);
197: DMCreateDS(dmProj);
199: DMClone(dma, &dmaProj);
200: PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, -1, &fe);
201: DMSetField(dmaProj, 0, NULL, (PetscObject)fe);
202: PetscFEDestroy(&fe);
203: DMCreateDS(dmaProj);
205: DMGetGlobalVector(dmProj, &u);
206: DMGetGlobalVector(dmaProj, &ua);
207: DMGetGlobalVector(dmaProj, &rhs);
208: DMGetGlobalVector(dmaProj, &uproj);
210: // Interpolate onto original mesh using dual basis
211: DMProjectFunction(dmProj, 0.0, funcs, NULL, INSERT_VALUES, u);
212: PetscObjectSetName((PetscObject)u, "Original");
213: VecViewFromOptions(u, NULL, "-orig_vec_view");
214: DMComputeL2Diff(dmProj, 0.0, funcs, NULL, u, &error);
215: PetscPrintf(PETSC_COMM_WORLD, "Original L2 Error: %g\n", (double)error);
216: // Interpolate onto NEW mesh using dual basis
217: DMProjectFunction(dmaProj, 0.0, funcs, NULL, INSERT_VALUES, ua);
218: PetscObjectSetName((PetscObject)ua, "Adapted");
219: VecViewFromOptions(ua, NULL, "-adapt_vec_view");
220: DMComputeL2Diff(dmaProj, 0.0, funcs, NULL, ua, &error);
221: PetscPrintf(PETSC_COMM_WORLD, "Adapted L2 Error: %g\n", (double)error);
222: // Interpolate between meshes using interpolation matrix
223: DMCreateInterpolation(dmProj, dmaProj, &Interp, &scaling);
224: MatInterpolate(Interp, u, ua);
225: MatDestroy(&Interp);
226: VecDestroy(&scaling);
227: PetscObjectSetName((PetscObject)ua, "Interpolation");
228: VecViewFromOptions(ua, NULL, "-interp_vec_view");
229: DMComputeL2Diff(dmaProj, 0.0, funcs, NULL, ua, &error);
230: PetscPrintf(PETSC_COMM_WORLD, "Interpolated L2 Error: %g\n", (double)error);
231: // L2 projection
232: DMCreateMassMatrix(dmaProj, dmaProj, &mass);
233: MatViewFromOptions(mass, NULL, "-mass_mat_view");
234: KSPCreate(PETSC_COMM_WORLD, &ksp);
235: KSPSetOperators(ksp, mass, mass);
236: KSPSetFromOptions(ksp);
237: // Compute rhs as M f, could also directly project the analytic function but we might not have it
238: DMCreateMassMatrix(dmProj, dmaProj, &mass2);
239: MatMult(mass2, u, rhs);
240: MatDestroy(&mass2);
241: KSPSolve(ksp, rhs, uproj);
242: PetscObjectSetName((PetscObject)uproj, "L_2 Projection");
243: VecViewFromOptions(uproj, NULL, "-proj_vec_view");
244: DMComputeL2Diff(dmaProj, 0.0, funcs, NULL, uproj, &error);
245: PetscPrintf(PETSC_COMM_WORLD, "Projected L2 Error: %g\n", (double)error);
246: KSPDestroy(&ksp);
247: MatDestroy(&mass);
248: DMRestoreGlobalVector(dmProj, &u);
249: DMRestoreGlobalVector(dmaProj, &ua);
250: DMRestoreGlobalVector(dmaProj, &rhs);
251: DMRestoreGlobalVector(dmaProj, &uproj);
252: DMDestroy(&dmProj);
253: DMDestroy(&dmaProj);
254: return 0;
255: }
257: int main(int argc, char *argv[])
258: {
259: DM dm;
260: AppCtx user; /* user-defined work context */
261: MPI_Comm comm;
262: DM dma, odm;
263: Vec metric;
264: PetscInt r;
267: PetscInitialize(&argc, &argv, NULL, help);
268: comm = PETSC_COMM_WORLD;
269: ProcessOptions(comm, &user);
270: CreateMesh(comm, &dm);
272: odm = dm;
273: DMPlexDistributeOverlap(odm, 1, NULL, &dm);
274: if (!dm) {
275: dm = odm;
276: } else DMDestroy(&odm);
278: for (r = 0; r < user.Nr; ++r) {
279: DMLabel label;
281: ComputeMetric(dm, &user, &metric);
282: DMGetLabel(dm, "marker", &label);
283: DMAdaptMetric(dm, metric, label, NULL, &dma);
284: VecDestroy(&metric);
285: PetscObjectSetName((PetscObject)dma, "DMadapt");
286: PetscObjectSetOptionsPrefix((PetscObject)dma, "adapt_");
287: DMViewFromOptions(dma, NULL, "-dm_view");
288: if (user.doL2) TestL2Projection(dm, dma, &user);
289: DMDestroy(&dm);
290: dm = dma;
291: }
292: PetscObjectSetOptionsPrefix((PetscObject)dm, "final_");
293: DMViewFromOptions(dm, NULL, "-dm_view");
294: DMDestroy(&dm);
295: PetscFinalize();
296: return 0;
297: }
299: /*TEST
301: build:
302: requires: pragmatic
304: testset:
305: args: -dm_plex_box_faces 4,4,4 -dm_adaptor pragmatic -met 2 -init_dm_view -adapt_dm_view -dm_adaptor pragmatic
307: test:
308: suffix: 2d
309: args: -dm_plex_separate_marker 0
310: test:
311: suffix: 2d_sep
312: args: -dm_plex_separate_marker 1
313: test:
314: suffix: 3d
315: args: -dm_plex_dim 3
317: # Pragmatic hangs for simple partitioner
318: testset:
319: requires: parmetis
320: args: -dm_plex_box_faces 2,2 -petscpartitioner_type parmetis -met 2 -init_dm_view -adapt_dm_view -dm_adaptor pragmatic
322: test:
323: suffix: 2d_parmetis_np2
324: nsize: 2
325: test:
326: suffix: 2d_parmetis_np4
327: nsize: 4
329: test:
330: requires: parmetis
331: suffix: 3d_parmetis_met0
332: nsize: 2
333: args: -dm_plex_dim 3 -dm_plex_box_faces 9,9,9 -dm_adaptor pragmatic -petscpartitioner_type parmetis \
334: -met 0 -hmin 0.01 -hmax 0.03 -init_dm_view -adapt_dm_view -dm_adaptor pragmatic
335: test:
336: requires: parmetis
337: suffix: 3d_parmetis_met2
338: nsize: 2
339: args: -dm_plex_box_faces 19,19 -dm_adaptor pragmatic -petscpartitioner_type parmetis \
340: -met 2 -hmax 0.5 -hmin 0.001 -init_dm_view -adapt_dm_view -dm_adaptor pragmatic
341: test:
342: suffix: proj2
343: args: -dm_plex_box_faces 2,2 -dm_plex_hash_location -dm_adaptor pragmatic -init_dm_view -adapt_dm_view -do_L2 \
344: -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu -dm_adaptor pragmatic
345: test:
346: suffix: proj4
347: args: -dm_plex_box_faces 4,4 -dm_plex_hash_location -dm_adaptor pragmatic -init_dm_view -adapt_dm_view -do_L2 \
348: -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu -dm_adaptor pragmatic
350: test:
351: suffix: 2d_met3
352: args: -dm_plex_box_faces 9,9 -met 3 -dm_adaptor pragmatic -init_dm_view -adapt_dm_view \
353: -dm_plex_metric_h_min 1.e-10 -dm_plex_metric_h_max 1.0e-01 -dm_plex_metric_a_max 1.0e+05 -dm_plex_metric_p 1.0 \
354: -dm_plex_metric_target_complexity 10000.0 -dm_adaptor pragmatic
356: TEST*/