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