Actual source code: ex4.c
1: #include <petscsnes.h>
2: #include <petscdmda.h>
4: static const char help[] = "Minimum surface area problem in 2D using DMDA.\n\
5: It solves an unconstrained minimization problem. This example is based on a \n\
6: problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and \n\
7: boundary values along the edges of the domain, the objective is to find the\n\
8: surface with the minimal area that satisfies the boundary conditions.\n\
9: \n\
10: The command line options are:\n\
11: -da_grid_x <nx>, where <nx> = number of grid points in the 1st coordinate direction\n\
12: -da_grid_y <ny>, where <ny> = number of grid points in the 2nd coordinate direction\n\
13: \n";
15: /*
16: User-defined application context - contains data needed by the
17: application-provided call-back routines, FormJacobianLocal() and
18: FormFunctionLocal().
19: */
21: typedef enum {
22: PROBLEM_ENNEPER,
23: PROBLEM_SINS,
24: } ProblemType;
25: static const char *const ProblemTypes[] = {"ENNEPER", "SINS", "ProblemType", "PROBLEM_", 0};
27: typedef struct {
28: PetscScalar *bottom, *top, *left, *right;
29: } AppCtx;
31: /* -------- User-defined Routines --------- */
33: static PetscErrorCode FormBoundaryConditions_Enneper(SNES, AppCtx **);
34: static PetscErrorCode FormBoundaryConditions_Sins(SNES, AppCtx **);
35: static PetscErrorCode DestroyBoundaryConditions(AppCtx **);
36: static PetscErrorCode FormObjectiveLocal(DMDALocalInfo *, PetscScalar **, PetscReal *, void *);
37: static PetscErrorCode FormFunctionLocal(DMDALocalInfo *, PetscScalar **, PetscScalar **, void *);
38: static PetscErrorCode FormJacobianLocal(DMDALocalInfo *, PetscScalar **, Mat, Mat, void *);
40: int main(int argc, char **argv)
41: {
42: Vec x;
43: SNES snes;
44: DM da;
45: ProblemType ptype = PROBLEM_ENNEPER;
46: PetscBool use_obj = PETSC_TRUE;
47: PetscReal bbox[4] = {0.};
48: AppCtx *user;
49: PetscErrorCode (*form_bc)(SNES, AppCtx **) = NULL;
51: PetscFunctionBeginUser;
52: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
54: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Minimal surface options", __FILE__);
55: PetscCall(PetscOptionsEnum("-problem_type", "Problem type", NULL, ProblemTypes, (PetscEnum)ptype, (PetscEnum *)&ptype, NULL));
56: PetscCall(PetscOptionsBool("-use_objective", "Use objective function", NULL, use_obj, &use_obj, NULL));
57: PetscOptionsEnd();
58: switch (ptype) {
59: case PROBLEM_ENNEPER:
60: bbox[0] = -0.5;
61: bbox[1] = 0.5;
62: bbox[2] = -0.5;
63: bbox[3] = 0.5;
64: form_bc = FormBoundaryConditions_Enneper;
65: break;
66: case PROBLEM_SINS:
67: bbox[0] = 0.0;
68: bbox[1] = 1.0;
69: bbox[2] = 0.0;
70: bbox[3] = 1.0;
71: form_bc = FormBoundaryConditions_Sins;
72: break;
73: }
75: /* Create distributed array to manage the 2d grid */
76: PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da));
77: PetscCall(DMSetFromOptions(da));
78: PetscCall(DMSetUp(da));
79: PetscCall(DMDASetUniformCoordinates(da, bbox[0], bbox[1], bbox[2], bbox[3], PETSC_DECIDE, PETSC_DECIDE));
81: /* Extract global vectors from DMDA; */
82: PetscCall(DMCreateGlobalVector(da, &x));
84: /* Create nonlinear solver context */
85: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
86: PetscCall(SNESSetDM(snes, da));
87: PetscCall((*form_bc)(snes, &user));
88: PetscCall(SNESSetApplicationContext(snes, user));
90: /* Set local callbacks */
91: if (use_obj) PetscCall(DMDASNESSetObjectiveLocal(da, (DMDASNESObjectiveFn *)FormObjectiveLocal, NULL));
92: PetscCall(DMDASNESSetFunctionLocal(da, INSERT_VALUES, (DMDASNESFunctionFn *)FormFunctionLocal, NULL));
93: PetscCall(DMDASNESSetJacobianLocal(da, (DMDASNESJacobianFn *)FormJacobianLocal, NULL));
95: /* Customize from command line */
96: PetscCall(SNESSetFromOptions(snes));
98: /* Solve the application */
99: PetscCall(SNESSolve(snes, NULL, x));
101: /* Free user-created data structures */
102: PetscCall(VecDestroy(&x));
103: PetscCall(SNESDestroy(&snes));
104: PetscCall(DMDestroy(&da));
105: PetscCall(DestroyBoundaryConditions(&user));
107: PetscCall(PetscFinalize());
108: return 0;
109: }
111: /* Compute objective function over the locally owned part of the mesh */
112: PetscErrorCode FormObjectiveLocal(DMDALocalInfo *info, PetscScalar **x, PetscReal *v, void *ptr)
113: {
114: AppCtx *user = (AppCtx *)ptr;
115: PetscInt mx = info->mx, my = info->my;
116: PetscInt xs = info->xs, xm = info->xm, ys = info->ys, ym = info->ym;
117: PetscInt i, j;
118: PetscScalar hx, hy;
119: PetscScalar f2, f4, d1, d2, d3, d4, xc, xl, xr, xt, xb;
120: PetscReal ft = 0, area;
122: PetscFunctionBeginUser;
123: PetscCheck(user, PetscObjectComm((PetscObject)info->da), PETSC_ERR_PLIB, "Missing application context");
124: hx = 1.0 / (mx + 1);
125: hy = 1.0 / (my + 1);
126: area = 0.5 * hx * hy;
127: for (j = ys; j < ys + ym; j++) {
128: for (i = xs; i < xs + xm; i++) {
129: xc = x[j][i];
130: xl = xr = xb = xt = xc;
132: if (i == 0) { /* left side */
133: xl = user->left[j + 1];
134: } else xl = x[j][i - 1];
136: if (j == 0) { /* bottom side */
137: xb = user->bottom[i + 1];
138: } else xb = x[j - 1][i];
140: if (i + 1 == mx) { /* right side */
141: xr = user->right[j + 1];
142: } else xr = x[j][i + 1];
144: if (j + 1 == 0 + my) { /* top side */
145: xt = user->top[i + 1];
146: } else xt = x[j + 1][i];
148: d1 = (xc - xl);
149: d2 = (xc - xr);
150: d3 = (xc - xt);
151: d4 = (xc - xb);
153: d1 /= hx;
154: d2 /= hx;
155: d3 /= hy;
156: d4 /= hy;
158: f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
159: f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
161: ft += PetscRealPart(f2 + f4);
162: }
163: }
165: /* Compute triangular areas along the border of the domain. */
166: if (xs == 0) { /* left side */
167: for (j = ys; j < ys + ym; j++) {
168: d3 = (user->left[j + 1] - user->left[j + 2]) / hy;
169: d2 = (user->left[j + 1] - x[j][0]) / hx;
170: ft += PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
171: }
172: }
173: if (ys == 0) { /* bottom side */
174: for (i = xs; i < xs + xm; i++) {
175: d2 = (user->bottom[i + 1] - user->bottom[i + 2]) / hx;
176: d3 = (user->bottom[i + 1] - x[0][i]) / hy;
177: ft += PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
178: }
179: }
180: if (xs + xm == mx) { /* right side */
181: for (j = ys; j < ys + ym; j++) {
182: d1 = (x[j][mx - 1] - user->right[j + 1]) / hx;
183: d4 = (user->right[j] - user->right[j + 1]) / hy;
184: ft += PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
185: }
186: }
187: if (ys + ym == my) { /* top side */
188: for (i = xs; i < xs + xm; i++) {
189: d1 = (x[my - 1][i] - user->top[i + 1]) / hy;
190: d4 = (user->top[i + 1] - user->top[i]) / hx;
191: ft += PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
192: }
193: }
194: if (ys == 0 && xs == 0) {
195: d1 = (user->left[0] - user->left[1]) / hy;
196: d2 = (user->bottom[0] - user->bottom[1]) / hx;
197: ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
198: }
199: if (ys + ym == my && xs + xm == mx) {
200: d1 = (user->right[ym + 1] - user->right[ym]) / hy;
201: d2 = (user->top[xm + 1] - user->top[xm]) / hx;
202: ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
203: }
204: ft *= area;
205: *v = ft;
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: /* Compute gradient over the locally owned part of the mesh */
210: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, PetscScalar **x, PetscScalar **g, void *ptr)
211: {
212: AppCtx *user = (AppCtx *)ptr;
213: PetscInt mx = info->mx, my = info->my;
214: PetscInt xs = info->xs, xm = info->xm, ys = info->ys, ym = info->ym;
215: PetscInt i, j;
216: PetscScalar hx, hy, hydhx, hxdhy;
217: PetscScalar f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
218: PetscScalar df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc;
220: PetscFunctionBeginUser;
221: PetscCheck(user, PetscObjectComm((PetscObject)info->da), PETSC_ERR_PLIB, "Missing application context");
222: hx = 1.0 / (mx + 1);
223: hy = 1.0 / (my + 1);
224: hydhx = hy / hx;
225: hxdhy = hx / hy;
227: for (j = ys; j < ys + ym; j++) {
228: for (i = xs; i < xs + xm; i++) {
229: xc = x[j][i];
230: xlt = xrb = xl = xr = xb = xt = xc;
232: if (i == 0) { /* left side */
233: xl = user->left[j + 1];
234: xlt = user->left[j + 2];
235: } else xl = x[j][i - 1];
237: if (j == 0) { /* bottom side */
238: xb = user->bottom[i + 1];
239: xrb = user->bottom[i + 2];
240: } else xb = x[j - 1][i];
242: if (i + 1 == mx) { /* right side */
243: xr = user->right[j + 1];
244: xrb = user->right[j];
245: } else xr = x[j][i + 1];
247: if (j + 1 == 0 + my) { /* top side */
248: xt = user->top[i + 1];
249: xlt = user->top[i];
250: } else xt = x[j + 1][i];
252: if (i > 0 && j + 1 < my) xlt = x[j + 1][i - 1]; /* left top side */
253: if (j > 0 && i + 1 < mx) xrb = x[j - 1][i + 1]; /* right bottom */
255: d1 = (xc - xl);
256: d2 = (xc - xr);
257: d3 = (xc - xt);
258: d4 = (xc - xb);
259: d5 = (xr - xrb);
260: d6 = (xrb - xb);
261: d7 = (xlt - xl);
262: d8 = (xt - xlt);
264: df1dxc = d1 * hydhx;
265: df2dxc = (d1 * hydhx + d4 * hxdhy);
266: df3dxc = d3 * hxdhy;
267: df4dxc = (d2 * hydhx + d3 * hxdhy);
268: df5dxc = d2 * hydhx;
269: df6dxc = d4 * hxdhy;
271: d1 /= hx;
272: d2 /= hx;
273: d3 /= hy;
274: d4 /= hy;
275: d5 /= hy;
276: d6 /= hx;
277: d7 /= hy;
278: d8 /= hx;
280: f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
281: f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
282: f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
283: f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
284: f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
285: f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);
287: df1dxc /= f1;
288: df2dxc /= f2;
289: df3dxc /= f3;
290: df4dxc /= f4;
291: df5dxc /= f5;
292: df6dxc /= f6;
294: g[j][i] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) / 2.0;
295: }
296: }
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: /* Compute Hessian over the locally owned part of the mesh */
301: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, PetscScalar **x, Mat H, Mat Hp, void *ptr)
302: {
303: AppCtx *user = (AppCtx *)ptr;
304: PetscInt mx = info->mx, my = info->my;
305: PetscInt xs = info->xs, xm = info->xm, ys = info->ys, ym = info->ym;
306: PetscInt i, j, k;
307: MatStencil row, col[7];
308: PetscScalar hx, hy, hydhx, hxdhy;
309: PetscScalar f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
310: PetscScalar hl, hr, ht, hb, hc, htl, hbr;
311: PetscScalar v[7];
313: PetscFunctionBeginUser;
314: PetscCheck(user, PetscObjectComm((PetscObject)info->da), PETSC_ERR_PLIB, "Missing application context");
315: hx = 1.0 / (mx + 1);
316: hy = 1.0 / (my + 1);
317: hydhx = hy / hx;
318: hxdhy = hx / hy;
320: for (j = ys; j < ys + ym; j++) {
321: for (i = xs; i < xs + xm; i++) {
322: xc = x[j][i];
323: xlt = xrb = xl = xr = xb = xt = xc;
325: /* Left */
326: if (i == 0) {
327: xl = user->left[j + 1];
328: xlt = user->left[j + 2];
329: } else xl = x[j][i - 1];
331: /* Bottom */
332: if (j == 0) {
333: xb = user->bottom[i + 1];
334: xrb = user->bottom[i + 2];
335: } else xb = x[j - 1][i];
337: /* Right */
338: if (i + 1 == mx) {
339: xr = user->right[j + 1];
340: xrb = user->right[j];
341: } else xr = x[j][i + 1];
343: /* Top */
344: if (j + 1 == my) {
345: xt = user->top[i + 1];
346: xlt = user->top[i];
347: } else xt = x[j + 1][i];
349: /* Top left */
350: if (i > 0 && j + 1 < my) xlt = x[j + 1][i - 1];
352: /* Bottom right */
353: if (j > 0 && i + 1 < mx) xrb = x[j - 1][i + 1];
355: d1 = (xc - xl) / hx;
356: d2 = (xc - xr) / hx;
357: d3 = (xc - xt) / hy;
358: d4 = (xc - xb) / hy;
359: d5 = (xrb - xr) / hy;
360: d6 = (xrb - xb) / hx;
361: d7 = (xlt - xl) / hy;
362: d8 = (xlt - xt) / hx;
364: f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
365: f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
366: f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
367: f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
368: f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
369: f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);
371: hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2);
372: hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4);
373: ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4);
374: hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2);
376: hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6);
377: htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3);
379: hc = hydhx * (1.0 + d7 * d7) / (f1 * f1 * f1) + hxdhy * (1.0 + d8 * d8) / (f3 * f3 * f3) + hydhx * (1.0 + d5 * d5) / (f5 * f5 * f5) + hxdhy * (1.0 + d6 * d6) / (f6 * f6 * f6) + (hxdhy * (1.0 + d1 * d1) + hydhx * (1.0 + d4 * d4) - 2.0 * d1 * d4) / (f2 * f2 * f2) + (hxdhy * (1.0 + d2 * d2) + hydhx * (1.0 + d3 * d3) - 2.0 * d2 * d3) / (f4 * f4 * f4);
381: hl /= 2.0;
382: hr /= 2.0;
383: ht /= 2.0;
384: hb /= 2.0;
385: hbr /= 2.0;
386: htl /= 2.0;
387: hc /= 2.0;
389: k = 0;
390: row.i = i;
391: row.j = j;
392: /* Bottom */
393: if (j > 0) {
394: v[k] = hb;
395: col[k].i = i;
396: col[k].j = j - 1;
397: k++;
398: }
400: /* Bottom right */
401: if (j > 0 && i < mx - 1) {
402: v[k] = hbr;
403: col[k].i = i + 1;
404: col[k].j = j - 1;
405: k++;
406: }
408: /* left */
409: if (i > 0) {
410: v[k] = hl;
411: col[k].i = i - 1;
412: col[k].j = j;
413: k++;
414: }
416: /* Centre */
417: v[k] = hc;
418: col[k].i = row.i;
419: col[k].j = row.j;
420: k++;
422: /* Right */
423: if (i < mx - 1) {
424: v[k] = hr;
425: col[k].i = i + 1;
426: col[k].j = j;
427: k++;
428: }
430: /* Top left */
431: if (i > 0 && j < my - 1) {
432: v[k] = htl;
433: col[k].i = i - 1;
434: col[k].j = j + 1;
435: k++;
436: }
438: /* Top */
439: if (j < my - 1) {
440: v[k] = ht;
441: col[k].i = i;
442: col[k].j = j + 1;
443: k++;
444: }
446: PetscCall(MatSetValuesStencil(Hp, 1, &row, k, col, v, INSERT_VALUES));
447: }
448: }
450: /* Assemble the matrix */
451: PetscCall(MatAssemblyBegin(Hp, MAT_FINAL_ASSEMBLY));
452: PetscCall(MatAssemblyEnd(Hp, MAT_FINAL_ASSEMBLY));
453: PetscFunctionReturn(PETSC_SUCCESS);
454: }
456: PetscErrorCode FormBoundaryConditions_Enneper(SNES snes, AppCtx **ouser)
457: {
458: PetscInt i, j, k, limit = 0, maxits = 5;
459: PetscInt mx, my;
460: PetscInt bsize = 0, lsize = 0, tsize = 0, rsize = 0;
461: PetscScalar one = 1.0, two = 2.0, three = 3.0;
462: PetscScalar det, hx, hy, xt = 0, yt = 0;
463: PetscReal fnorm, tol = 1e-10;
464: PetscScalar u1, u2, nf1, nf2, njac11, njac12, njac21, njac22;
465: PetscScalar b = -0.5, t = 0.5, l = -0.5, r = 0.5;
466: PetscScalar *boundary;
467: AppCtx *user;
468: DM da;
470: PetscFunctionBeginUser;
471: PetscCall(SNESGetDM(snes, &da));
472: PetscCall(PetscNew(&user));
473: *ouser = user;
474: PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &mx, &my, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
475: bsize = mx + 2;
476: lsize = my + 2;
477: rsize = my + 2;
478: tsize = mx + 2;
480: PetscCall(PetscMalloc1(bsize, &user->bottom));
481: PetscCall(PetscMalloc1(tsize, &user->top));
482: PetscCall(PetscMalloc1(lsize, &user->left));
483: PetscCall(PetscMalloc1(rsize, &user->right));
485: hx = 1.0 / (mx + 1.0);
486: hy = 1.0 / (my + 1.0);
488: for (j = 0; j < 4; j++) {
489: if (j == 0) {
490: yt = b;
491: xt = l;
492: limit = bsize;
493: boundary = user->bottom;
494: } else if (j == 1) {
495: yt = t;
496: xt = l;
497: limit = tsize;
498: boundary = user->top;
499: } else if (j == 2) {
500: yt = b;
501: xt = l;
502: limit = lsize;
503: boundary = user->left;
504: } else { /* if (j==3) */
505: yt = b;
506: xt = r;
507: limit = rsize;
508: boundary = user->right;
509: }
511: for (i = 0; i < limit; i++) {
512: u1 = xt;
513: u2 = -yt;
514: for (k = 0; k < maxits; k++) {
515: nf1 = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt;
516: nf2 = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt;
517: fnorm = PetscRealPart(PetscSqrtScalar(nf1 * nf1 + nf2 * nf2));
518: if (fnorm <= tol) break;
519: njac11 = one + u2 * u2 - u1 * u1;
520: njac12 = two * u1 * u2;
521: njac21 = -two * u1 * u2;
522: njac22 = -one - u1 * u1 + u2 * u2;
523: det = njac11 * njac22 - njac21 * njac12;
524: u1 = u1 - (njac22 * nf1 - njac12 * nf2) / det;
525: u2 = u2 - (njac11 * nf2 - njac21 * nf1) / det;
526: }
528: boundary[i] = u1 * u1 - u2 * u2;
529: if (j == 0 || j == 1) xt = xt + hx;
530: else yt = yt + hy; /* if (j==2 || j==3) */
531: }
532: }
533: PetscFunctionReturn(PETSC_SUCCESS);
534: }
536: PetscErrorCode DestroyBoundaryConditions(AppCtx **ouser)
537: {
538: AppCtx *user = *ouser;
540: PetscFunctionBeginUser;
541: PetscCall(PetscFree(user->bottom));
542: PetscCall(PetscFree(user->top));
543: PetscCall(PetscFree(user->left));
544: PetscCall(PetscFree(user->right));
545: PetscCall(PetscFree(*ouser));
546: PetscFunctionReturn(PETSC_SUCCESS);
547: }
549: PetscErrorCode FormBoundaryConditions_Sins(SNES snes, AppCtx **ouser)
550: {
551: PetscInt i, j;
552: PetscInt mx, my;
553: PetscInt limit, bsize = 0, lsize = 0, tsize = 0, rsize = 0;
554: PetscScalar hx, hy, xt = 0, yt = 0;
555: PetscScalar b = 0.0, t = 1.0, l = 0.0, r = 1.0;
556: PetscScalar *boundary;
557: AppCtx *user;
558: DM da;
559: PetscReal pi2 = 2 * PETSC_PI;
561: PetscFunctionBeginUser;
562: PetscCall(SNESGetDM(snes, &da));
563: PetscCall(PetscNew(&user));
564: *ouser = user;
565: PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &mx, &my, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
566: bsize = mx + 2;
567: lsize = my + 2;
568: rsize = my + 2;
569: tsize = mx + 2;
571: PetscCall(PetscMalloc1(bsize, &user->bottom));
572: PetscCall(PetscMalloc1(tsize, &user->top));
573: PetscCall(PetscMalloc1(lsize, &user->left));
574: PetscCall(PetscMalloc1(rsize, &user->right));
576: hx = 1.0 / (mx + 1.0);
577: hy = 1.0 / (my + 1.0);
579: for (j = 0; j < 4; j++) {
580: if (j == 0) {
581: yt = b;
582: xt = l;
583: limit = bsize;
584: boundary = user->bottom;
585: } else if (j == 1) {
586: yt = t;
587: xt = l;
588: limit = tsize;
589: boundary = user->top;
590: } else if (j == 2) {
591: yt = b;
592: xt = l;
593: limit = lsize;
594: boundary = user->left;
595: } else { /* if (j==3) */
596: yt = b;
597: xt = r;
598: limit = rsize;
599: boundary = user->right;
600: }
602: for (i = 0; i < limit; i++) {
603: if (j == 0) { /* bottom */
604: boundary[i] = -0.5 * PetscSinReal(pi2 * xt);
605: } else if (j == 1) { /* top */
606: boundary[i] = 0.5 * PetscSinReal(pi2 * xt);
607: } else if (j == 2) { /* left */
608: boundary[i] = -0.5 * PetscSinReal(pi2 * yt);
609: } else { /* right */
610: boundary[i] = 0.5 * PetscSinReal(pi2 * yt);
611: }
612: if (j == 0 || j == 1) xt = xt + hx;
613: else yt = yt + hy;
614: }
615: }
616: PetscFunctionReturn(PETSC_SUCCESS);
617: }
619: /*TEST
621: build:
622: requires: !complex
624: test:
625: requires: !single
626: filter: sed -e "s/CONVERGED_FNORM_ABS/CONVERGED_FNORM_RELATIVE/g"
627: suffix: qn_nasm
628: args: -snes_type qn -snes_npc_side {{left right}separate output} -npc_snes_type nasm -snes_converged_reason -da_local_subdomains 4 -use_objective {{0 1}separate output}
630: TEST*/