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