1 static char help[] = "Tests mesh adaptation with DMPlex and pragmatic.\n";
2
3 #include <petsc/private/dmpleximpl.h>
4
5 #include <petscsnes.h>
6
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;
13
14 /*
15 Classic hyperbolic sensor function for testing multi-scale anisotropic mesh adaptation:
16
17 f:[-1, 1]x[-1, 1] \to R,
18 f(x, y) = sin(50xy)/100 if |xy| > 2\pi/50 else sin(50xy)
19
20 (mapped to have domain [0,1] x [0,1] in this case).
21 */
sensor(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar u[],PetscCtx ctx)22 static PetscErrorCode sensor(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], PetscCtx ctx)
23 {
24 const PetscReal xref = 2. * x[0] - 1.;
25 const PetscReal yref = 2. * x[1] - 1.;
26 const PetscReal xy = xref * yref;
27
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 }
33
ProcessOptions(MPI_Comm comm,AppCtx * options)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;
42
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 }
52
CreateMesh(MPI_Comm comm,DM * dm)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 }
63
ComputeMetricSensor(DM dm,AppCtx * user,Vec * metric)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;
72
73 PetscFunctionBegin;
74 PetscCall(DMGetDimension(dm, &dim));
75 PetscCall(DMPlexIsSimplex(dm, &simplex));
76
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"));
85
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"));
96
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"));
103
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 }
116
ComputeMetric(DM dm,AppCtx * user,Vec * metric)117 static PetscErrorCode ComputeMetric(DM dm, AppCtx *user, Vec *metric)
118 {
119 PetscReal lambda = 1 / (user->hmax * user->hmax);
120
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;
134
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;
145
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 }
172
linear(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)173 static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
174 {
175 u[0] = x[0] + x[1];
176 return 0;
177 }
178
TestL2Projection(DM dm,DM dma,AppCtx * user)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;
190
191 PetscFunctionBeginUser;
192 PetscCall(DMGetDimension(dm, &dim));
193 PetscCall(DMPlexIsSimplex(dm, &simplex));
194
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));
200
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));
206
207 PetscCall(DMGetGlobalVector(dmProj, &u));
208 PetscCall(DMGetGlobalVector(dmaProj, &ua));
209 PetscCall(DMGetGlobalVector(dmaProj, &rhs));
210 PetscCall(DMGetGlobalVector(dmaProj, &uproj));
211
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 }
258
main(int argc,char * argv[])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;
267
268 PetscFunctionBeginUser;
269 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
270 comm = PETSC_COMM_WORLD;
271 PetscCall(ProcessOptions(comm, &user));
272 PetscCall(CreateMesh(comm, &dm));
273
274 odm = dm;
275 PetscCall(DMPlexDistributeOverlap(odm, 1, NULL, &dm));
276 if (!dm) {
277 dm = odm;
278 } else PetscCall(DMDestroy(&odm));
279
280 for (r = 0; r < user.Nr; ++r) {
281 DMLabel label;
282
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 }
300
301 /*TEST
302
303 build:
304 requires: pragmatic
305
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
308
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
318
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
323
324 test:
325 suffix: 2d_parmetis_np2
326 nsize: 2
327 test:
328 suffix: 2d_parmetis_np4
329 nsize: 4
330
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
351
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
357
358 TEST*/
359