Actual source code: ex16.c
petsc-3.3-p6 2013-02-11
1: #include <petscsnes.h>
2: #include <petscdmda.h>
4: static char help[]=
5: "This example is an implementation of minimal surface area with \n\
6: a plate problem from the TAO package (examples/plate2.c) \n\
7: This example is based on a problem from the MINPACK-2 test suite.\n\
8: Given a rectangular 2-D domain, boundary values along the edges of \n\
9: the domain, and a plate represented by lower boundary conditions, \n\
10: the objective is to find the surface with the minimal area that \n\
11: satisfies the boundary conditions.\n\
12: The command line options are:\n\
13: -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction\n\
14: -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction\n\
15: -bheight <ht>, where <ht> = height of the plate\n\
16: -start <st>, where <st> =0 for zero vector, <st> != 0 \n\
17: for an average of the boundary conditions\n\n";
19: /*
20: User-defined application context - contains data needed by the
21: application-provided call-back routines, FormJacobian() and
22: FormFunction().
23: */
25: typedef struct {
26: DM da;
27: Vec Bottom, Top, Left, Right;
28: PetscScalar bheight;
29: PetscInt mx,my,bmx,bmy;
30: } AppCtx;
33: /* -------- User-defined Routines --------- */
35: PetscErrorCode MSA_BoundaryConditions(AppCtx *);
36: PetscErrorCode MSA_InitialPoint(AppCtx *, Vec);
37: PetscErrorCode MSA_Plate(Vec,Vec,void*);
38: PetscErrorCode FormGradient(SNES, Vec, Vec, void *);
39: PetscErrorCode FormJacobian(SNES, Vec, Mat *, Mat*, MatStructure*,void *);
43: int main(int argc, char **argv)
44: {
45: PetscErrorCode info; /* used to check for functions returning nonzeros */
46: Vec x,r; /* solution and residual vectors */
47: Vec xl,xu; /* Bounds on the variables */
48: SNES snes; /* nonlinear solver context */
49: Mat J; /* Jacobian matrix */
50: AppCtx user; /* user-defined work context */
51: PetscBool flg;
53: /* Initialize PETSc */
54: PetscInitialize(&argc, &argv, (char *)0, help );
56: #if defined(PETSC_USE_COMPLEX)
57: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This example does not work for scalar type complex\n");
58: #endif
60: /* Create distributed array to manage the 2d grid */
61: info = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,-10,-10,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,&user.da);CHKERRQ(info);
62: info = DMDAGetInfo(user.da,PETSC_IGNORE,&user.mx,&user.my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(info);
64: user.bheight=0.1;
65: info = PetscOptionsGetScalar(PETSC_NULL,"-bheight",&user.bheight,&flg); CHKERRQ(info);
67: user.bmx = user.mx/2; user.bmy = user.my/2;
68: info = PetscOptionsGetInt(PETSC_NULL,"-bmx",&user.bmx,&flg); CHKERRQ(info);
69: info = PetscOptionsGetInt(PETSC_NULL,"-bmy",&user.bmy,&flg); CHKERRQ(info);
71: PetscPrintf(PETSC_COMM_WORLD,"\n---- Minimum Surface Area With Plate Problem -----\n");
72: PetscPrintf(PETSC_COMM_WORLD,"mx:%d, my:%d, bmx:%d, bmy:%d, height:%4.2f\n",
73: user.mx,user.my,user.bmx,user.bmy,user.bheight);
75: /* Extract global vectors from DMDA; */
76: info = DMCreateGlobalVector(user.da,&x);CHKERRQ(info);
77: info = VecDuplicate(x, &r); CHKERRQ(info);
79: info = DMCreateMatrix(user.da,MATAIJ,&J);CHKERRQ(info);
81: /* Create nonlinear solver context */
82: info = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(info);
84: /* Set function evaluation and Jacobian evaluation routines */
85: info = SNESSetDM(snes,user.da);CHKERRQ(info);
86: info = SNESSetFunction(snes,r,FormGradient,&user);CHKERRQ(info);
87: info = SNESSetJacobian(snes,J,J,FormJacobian,&user);CHKERRQ(info);
89: /* Set the boundary conditions */
90: info = MSA_BoundaryConditions(&user); CHKERRQ(info);
92: /* Set initial solution guess */
93: info = MSA_InitialPoint(&user, x); CHKERRQ(info);
95: info = SNESSetFromOptions(snes);CHKERRQ(info);
97: /* Set Bounds on variables */
98: info = VecDuplicate(x, &xl); CHKERRQ(info);
99: info = VecDuplicate(x, &xu); CHKERRQ(info);
100: info = MSA_Plate(xl,xu,&user);CHKERRQ(info);
102: info = SNESVISetVariableBounds(snes,xl,xu);CHKERRQ(info);
104: /* Solve the application */
105: info = SNESSolve(snes,PETSC_NULL,x);CHKERRQ(info);
107: info = PetscOptionsHasName(PETSC_NULL,"-view_sol",&flg);CHKERRQ(info);
108: if (flg) { info = VecView(x,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info); }
110: /* Free memory */
111: info = VecDestroy(&x); CHKERRQ(info);
112: info = VecDestroy(&xl); CHKERRQ(info);
113: info = VecDestroy(&xu); CHKERRQ(info);
114: info = VecDestroy(&r); CHKERRQ(info);
115: info = MatDestroy(&J); CHKERRQ(info);
116: info = SNESDestroy(&snes); CHKERRQ(info);
118: /* Free user-created data structures */
119: info = DMDestroy(&user.da);CHKERRQ(info);
120: info = VecDestroy(&user.Bottom); CHKERRQ(info);
121: info = VecDestroy(&user.Top); CHKERRQ(info);
122: info = VecDestroy(&user.Left); CHKERRQ(info);
123: info = VecDestroy(&user.Right); CHKERRQ(info);
125: info = PetscFinalize();
127: return 0;
128: }
130: /* -------------------------------------------------------------------- */
134: /* FormGradient - Evaluates gradient of f.
136: Input Parameters:
137: . snes - the SNES context
138: . X - input vector
139: . ptr - optional user-defined context, as set by SNESSetFunction()
140:
141: Output Parameters:
142: . G - vector containing the newly evaluated gradient
143: */
144: PetscErrorCode FormGradient(SNES snes, Vec X, Vec G, void *ptr){
145: AppCtx *user = (AppCtx *) ptr;
146: int info;
147: PetscInt i,j;
148: PetscInt mx=user->mx, my=user->my;
149: PetscScalar hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
150: PetscScalar f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
151: PetscScalar df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
152: PetscScalar **g, **x;
153: PetscInt xs,xm,ys,ym;
154: Vec localX;
155: PetscScalar *top,*bottom,*left,*right;
158: /* Initialize vector to zero */
159: info = VecSet(G,0.0);CHKERRQ(info);
161: /* Get local vector */
162: info = DMGetLocalVector(user->da,&localX);CHKERRQ(info);
163: info = VecGetArray(user->Top,&top); CHKERRQ(info);
164: info = VecGetArray(user->Bottom,&bottom); CHKERRQ(info);
165: info = VecGetArray(user->Left,&left); CHKERRQ(info);
166: info = VecGetArray(user->Right,&right); CHKERRQ(info);
168: /* Get ghost points */
169: info = DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);CHKERRQ(info);
170: info = DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);CHKERRQ(info);
171: /* Get pointers to local vector data */
172: info = DMDAVecGetArray(user->da,localX, &x); CHKERRQ(info);
173: info = DMDAVecGetArray(user->da,G, &g); CHKERRQ(info);
175: info = DMDAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);CHKERRQ(info);
176: /* Compute function over the locally owned part of the mesh */
177: for (j=ys; j < ys+ym; j++){
178: for (i=xs; i< xs+xm; i++){
179:
180: xc = x[j][i];
181: xlt=xrb=xl=xr=xb=xt=xc;
182:
183: if (i==0){ /* left side */
184: xl= left[j-ys+1];
185: xlt = left[j-ys+2];
186: } else {
187: xl = x[j][i-1];
188: }
190: if (j==0){ /* bottom side */
191: xb=bottom[i-xs+1];
192: xrb = bottom[i-xs+2];
193: } else {
194: xb = x[j-1][i];
195: }
196:
197: if (i+1 == mx){ /* right side */
198: xr=right[j-ys+1];
199: xrb = right[j-ys];
200: } else {
201: xr = x[j][i+1];
202: }
204: if (j+1==0+my){ /* top side */
205: xt=top[i-xs+1];
206: xlt = top[i-xs];
207: }else {
208: xt = x[j+1][i];
209: }
211: if (i>0 && j+1<my){ /* left top side */
212: xlt = x[j+1][i-1];
213: }
214: if (j>0 && i+1<mx){ /* right bottom */
215: xrb = x[j-1][i+1];
216: }
218: d1 = (xc-xl);
219: d2 = (xc-xr);
220: d3 = (xc-xt);
221: d4 = (xc-xb);
222: d5 = (xr-xrb);
223: d6 = (xrb-xb);
224: d7 = (xlt-xl);
225: d8 = (xt-xlt);
226:
227: df1dxc = d1*hydhx;
228: df2dxc = ( d1*hydhx + d4*hxdhy );
229: df3dxc = d3*hxdhy;
230: df4dxc = ( d2*hydhx + d3*hxdhy );
231: df5dxc = d2*hydhx;
232: df6dxc = d4*hxdhy;
234: d1 /= hx;
235: d2 /= hx;
236: d3 /= hy;
237: d4 /= hy;
238: d5 /= hy;
239: d6 /= hx;
240: d7 /= hy;
241: d8 /= hx;
243: f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
244: f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
245: f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
246: f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
247: f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
248: f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);
249:
250: df1dxc /= f1;
251: df2dxc /= f2;
252: df3dxc /= f3;
253: df4dxc /= f4;
254: df5dxc /= f5;
255: df6dxc /= f6;
257: g[j][i] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc )/2.0;
258:
259: }
260: }
261:
262: /* Restore vectors */
263: info = DMDAVecRestoreArray(user->da,localX, &x); CHKERRQ(info);
264: info = DMDAVecRestoreArray(user->da,G, &g); CHKERRQ(info);
265: info = DMRestoreLocalVector(user->da,&localX);CHKERRQ(info);
267: info = VecRestoreArray(user->Left,&left); CHKERRQ(info);
268: info = VecRestoreArray(user->Top,&top); CHKERRQ(info);
269: info = VecRestoreArray(user->Bottom,&bottom); CHKERRQ(info);
270: info = VecRestoreArray(user->Right,&right); CHKERRQ(info);
272: info = PetscLogFlops(67*mx*my); CHKERRQ(info);
273: return(0);
274: }
276: /* ------------------------------------------------------------------- */
279: /*
280: FormJacobian - Evaluates Jacobian matrix.
282: Input Parameters:
283: . snes - SNES context
284: . X - input vector
285: . ptr - optional user-defined context, as set by SNESSetJacobian()
287: Output Parameters:
288: . tH - Jacobian matrix
290: */
291: PetscErrorCode FormJacobian(SNES snes, Vec X, Mat *tH, Mat* tHPre, MatStructure* flag, void *ptr)
292: {
293: AppCtx *user = (AppCtx *) ptr;
294: Mat H = *tH;
295: PetscErrorCode info;
296: PetscInt i,j,k;
297: PetscInt mx=user->mx, my=user->my;
298: MatStencil row,col[7];
299: PetscScalar hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
300: PetscScalar f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
301: PetscScalar hl,hr,ht,hb,hc,htl,hbr;
302: PetscScalar **x, v[7];
303: PetscBool assembled;
304: PetscInt xs,xm,ys,ym;
305: Vec localX;
306: PetscScalar *top,*bottom,*left,*right;
309: /* Set various matrix options */
310: info = MatAssembled(H,&assembled); CHKERRQ(info);
311: if (assembled){info = MatZeroEntries(H); CHKERRQ(info);}
312: *flag=SAME_NONZERO_PATTERN;
314: /* Get local vectors */
315: info = DMGetLocalVector(user->da,&localX);CHKERRQ(info);
316: info = VecGetArray(user->Top,&top); CHKERRQ(info);
317: info = VecGetArray(user->Bottom,&bottom); CHKERRQ(info);
318: info = VecGetArray(user->Left,&left); CHKERRQ(info);
319: info = VecGetArray(user->Right,&right); CHKERRQ(info);
321: /* Get ghost points */
322: info = DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);CHKERRQ(info);
323: info = DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);CHKERRQ(info);
324:
325: /* Get pointers to vector data */
326: info = DMDAVecGetArray(user->da,localX, &x); CHKERRQ(info);
328: info = DMDAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);CHKERRQ(info);
329: /* Compute Jacobian over the locally owned part of the mesh */
330: for (j=ys; j< ys+ym; j++){
331: for (i=xs; i< xs+xm; i++){
332: xc = x[j][i];
333: xlt=xrb=xl=xr=xb=xt=xc;
335: /* Left */
336: if (i==0){
337: xl= left[j+1];
338: xlt = left[j+2];
339: } else {
340: xl = x[j][i-1];
341: }
342:
343: /* Bottom */
344: if (j==0){
345: xb=bottom[i+1];
346: xrb = bottom[i+2];
347: } else {
348: xb = x[j-1][i];
349: }
350:
351: /* Right */
352: if (i+1 == mx){
353: xr=right[j+1];
354: xrb = right[j];
355: } else {
356: xr = x[j][i+1];
357: }
359: /* Top */
360: if (j+1==my){
361: xt=top[i+1];
362: xlt = top[i];
363: }else {
364: xt = x[j+1][i];
365: }
367: /* Top left */
368: if (i>0 && j+1<my){
369: xlt = x[j+1][i-1];
370: }
372: /* Bottom right */
373: if (j>0 && i+1<mx){
374: xrb = x[j-1][i+1];
375: }
377: d1 = (xc-xl)/hx;
378: d2 = (xc-xr)/hx;
379: d3 = (xc-xt)/hy;
380: d4 = (xc-xb)/hy;
381: d5 = (xrb-xr)/hy;
382: d6 = (xrb-xb)/hx;
383: d7 = (xlt-xl)/hy;
384: d8 = (xlt-xt)/hx;
385:
386: f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
387: f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
388: f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
389: f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
390: f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
391: f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);
394: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
395: (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
396: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
397: (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
398: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
399: (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
400: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
401: (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);
403: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
404: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);
406: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
407: hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
408: (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
409: (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);
411: hl/=2.0; hr/=2.0; ht/=2.0; hb/=2.0; hbr/=2.0; htl/=2.0; hc/=2.0;
413: k=0;
414: row.i = i;row.j= j;
415: /* Bottom */
416: if (j>0){
417: v[k]=hb;
418: col[k].i = i; col[k].j=j-1; k++;
419: }
420:
421: /* Bottom right */
422: if (j>0 && i < mx -1){
423: v[k]=hbr;
424: col[k].i = i+1; col[k].j = j-1; k++;
425: }
426:
427: /* left */
428: if (i>0){
429: v[k]= hl;
430: col[k].i = i-1; col[k].j = j; k++;
431: }
432:
433: /* Centre */
434: v[k]= hc; col[k].i= row.i; col[k].j = row.j; k++;
435:
436: /* Right */
437: if (i < mx-1 ){
438: v[k]= hr;
439: col[k].i= i+1; col[k].j = j;k++;
440: }
441:
442: /* Top left */
443: if (i>0 && j < my-1 ){
444: v[k]= htl;
445: col[k].i = i-1;col[k].j = j+1; k++;
446: }
447:
448: /* Top */
449: if (j < my-1 ){
450: v[k]= ht;
451: col[k].i = i; col[k].j = j+1; k++;
452: }
453:
454: info = MatSetValuesStencil(H,1,&row,k,col,v,INSERT_VALUES);
455: CHKERRQ(info);
456: }
457: }
459: info = VecRestoreArray(user->Left,&left); CHKERRQ(info);
460: info = VecRestoreArray(user->Top,&top); CHKERRQ(info);
461: info = VecRestoreArray(user->Bottom,&bottom); CHKERRQ(info);
462: info = VecRestoreArray(user->Right,&right); CHKERRQ(info);
464: /* Assemble the matrix */
465: info = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
466: info = DMDAVecRestoreArray(user->da,localX,&x);CHKERRQ(info);
467: info = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
468: info = DMRestoreLocalVector(user->da,&localX);CHKERRQ(info);
470: info = PetscLogFlops(199*mx*my); CHKERRQ(info);
471: return(0);
472: }
474: /* ------------------------------------------------------------------- */
477: /*
478: MSA_BoundaryConditions - Calculates the boundary conditions for
479: the region.
481: Input Parameter:
482: . user - user-defined application context
484: Output Parameter:
485: . user - user-defined application context
486: */
487: PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
488: {
489: PetscErrorCode info;
490: PetscInt i,j,k,limit=0,maxits=5;
491: PetscInt mx=user->mx,my=user->my;
492: PetscInt xs,ys,xm,ym;
493: PetscInt bsize=0, lsize=0, tsize=0, rsize=0;
494: PetscScalar one=1.0, two=2.0, three=3.0, tol=1e-10;
495: PetscScalar fnorm,det,hx,hy,xt=0,yt=0;
496: PetscScalar u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
497: PetscScalar b=-0.5, t=0.5, l=-0.5, r=0.5;
498: PetscScalar *boundary;
499: Vec Bottom,Top,Right,Left;
500: PetscScalar scl=1.0;
501: PetscBool flg;
505: info = DMDAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
507: bsize=xm+2; lsize=ym+2; rsize=ym+2; tsize=xm+2;
509: info = VecCreateMPI(PETSC_COMM_WORLD,bsize,PETSC_DECIDE,&Bottom); CHKERRQ(info);
510: info = VecCreateMPI(PETSC_COMM_WORLD,tsize,PETSC_DECIDE,&Top); CHKERRQ(info);
511: info = VecCreateMPI(PETSC_COMM_WORLD,lsize,PETSC_DECIDE,&Left); CHKERRQ(info);
512: info = VecCreateMPI(PETSC_COMM_WORLD,rsize,PETSC_DECIDE,&Right); CHKERRQ(info);
514: user->Top=Top;
515: user->Left=Left;
516: user->Bottom=Bottom;
517: user->Right=Right;
519: hx= (r-l)/(mx+1); hy=(t-b)/(my+1);
521: for (j=0; j<4; j++){
522: if (j==0){
523: yt=b;
524: xt=l+hx*xs;
525: limit=bsize;
526: info = VecGetArray(Bottom,&boundary);CHKERRQ(info);
527: } else if (j==1){
528: yt=t;
529: xt=l+hx*xs;
530: limit=tsize;
531: info = VecGetArray(Top,&boundary);CHKERRQ(info);
532: } else if (j==2){
533: yt=b+hy*ys;
534: xt=l;
535: limit=lsize;
536: info = VecGetArray(Left,&boundary); CHKERRQ(info);
537: } else { // if (j==3)
538: yt=b+hy*ys;
539: xt=r;
540: limit=rsize;
541: info = VecGetArray(Right,&boundary);CHKERRQ(info);
542: }
544: for (i=0; i<limit; i++){
545: u1=xt;
546: u2=-yt;
547: for (k=0; k<maxits; k++){
548: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
549: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
550: fnorm=PetscSqrtReal(nf1*nf1+nf2*nf2);
551: if (fnorm <= tol) break;
552: njac11=one+u2*u2-u1*u1;
553: njac12=two*u1*u2;
554: njac21=-two*u1*u2;
555: njac22=-one - u1*u1 + u2*u2;
556: det = njac11*njac22-njac21*njac12;
557: u1 = u1-(njac22*nf1-njac12*nf2)/det;
558: u2 = u2-(njac11*nf2-njac21*nf1)/det;
559: }
561: boundary[i]=u1*u1-u2*u2;
562: if (j==0 || j==1) {
563: xt=xt+hx;
564: } else { // if (j==2 || j==3)
565: yt=yt+hy;
566: }
567: }
569: if (j==0){
570: info = VecRestoreArray(Bottom,&boundary); CHKERRQ(info);
571: } else if (j==1){
572: info = VecRestoreArray(Top,&boundary); CHKERRQ(info);
573: } else if (j==2){
574: info = VecRestoreArray(Left,&boundary); CHKERRQ(info);
575: } else if (j==3){
576: info = VecRestoreArray(Right,&boundary); CHKERRQ(info);
577: }
579: }
581: /* Scale the boundary if desired */
583: info = PetscOptionsGetReal(PETSC_NULL,"-bottom",&scl,&flg);
584: CHKERRQ(info);
585: if (flg){
586: info = VecScale(Bottom, scl); CHKERRQ(info);
587: }
588:
589: info = PetscOptionsGetReal(PETSC_NULL,"-top",&scl,&flg);
590: CHKERRQ(info);
591: if (flg){
592: info = VecScale(Top, scl); CHKERRQ(info);
593: }
594:
595: info = PetscOptionsGetReal(PETSC_NULL,"-right",&scl,&flg);
596: CHKERRQ(info);
597: if (flg){
598: info = VecScale(Right, scl); CHKERRQ(info);
599: }
600:
601: info = PetscOptionsGetReal(PETSC_NULL,"-left",&scl,&flg);
602: CHKERRQ(info);
603: if (flg){
604: info = VecScale(Left, scl); CHKERRQ(info);
605: }
607: return(0);
608: }
610: /* ------------------------------------------------------------------- */
613: /*
614: MSA_InitialPoint - Calculates the initial guess in one of three ways.
616: Input Parameters:
617: . user - user-defined application context
618: . X - vector for initial guess
620: Output Parameters:
621: . X - newly computed initial guess
622: */
623: PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
624: {
625: PetscErrorCode info;
626: PetscInt start=-1,i,j;
627: PetscScalar zero=0.0;
628: PetscBool flg;
629: PetscScalar *left,*right,*bottom,*top;
632: info = PetscOptionsGetInt(PETSC_NULL,"-start",&start,&flg); CHKERRQ(info);
634: if (flg && start==0){ /* The zero vector is reasonable */
635:
636: info = VecSet(X, zero); CHKERRQ(info);
637: /* PLogInfo(user,"Min. Surface Area Problem: Start with 0 vector \n"); */
640: } else { /* Take an average of the boundary conditions */
641: PetscInt mx=user->mx,my=user->my;
642: PetscScalar **x;
643: PetscInt xs,xm,ys,ym;
644:
645: info = VecGetArray(user->Top,&top); CHKERRQ(info);
646: info = VecGetArray(user->Bottom,&bottom); CHKERRQ(info);
647: info = VecGetArray(user->Left,&left); CHKERRQ(info);
648: info = VecGetArray(user->Right,&right); CHKERRQ(info);
650: /* Get pointers to vector data */
651: info = DMDAVecGetArray(user->da,X,&x); CHKERRQ(info);
652: info = DMDAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);CHKERRQ(info);
654: /* Perform local computations */
655: for (j=ys; j<ys+ym; j++){
656: for (i=xs; i< xs+xm; i++){
657: x[j][i] = ( (j+1)*bottom[i-xs+1]/my+(my-j+1)*top[i-xs+1]/(my+2)+
658: (i+1)*left[j-ys+1]/mx+(mx-i+1)*right[j-ys+1]/(mx+2))/2.0;
659: }
660: }
661:
662: /* Restore vectors */
663: info = DMDAVecRestoreArray(user->da,X,&x); CHKERRQ(info);
664: info = VecRestoreArray(user->Left,&left); CHKERRQ(info);
665: info = VecRestoreArray(user->Top,&top); CHKERRQ(info);
666: info = VecRestoreArray(user->Bottom,&bottom); CHKERRQ(info);
667: info = VecRestoreArray(user->Right,&right); CHKERRQ(info);
669: }
670: return(0);
671: }
675: /*
676: MSA_Plate - Calculates an obstacle for surface to stretch over.
677: */
678: PetscErrorCode MSA_Plate(Vec XL,Vec XU,void *ctx)
679: {
680: AppCtx *user=(AppCtx *)ctx;
681: PetscErrorCode info;
682: PetscInt i,j;
683: PetscInt xs,ys,xm,ym;
684: PetscInt mx=user->mx, my=user->my, bmy, bmx;
685: PetscScalar t1,t2,t3;
686: PetscScalar **xl;
687: PetscScalar lb=-SNES_VI_INF, ub=SNES_VI_INF;
688: PetscBool cylinder;
690: user->bmy = PetscMax(0,user->bmy);user->bmy = PetscMin(my,user->bmy);
691: user->bmx = PetscMax(0,user->bmx);user->bmx = PetscMin(mx,user->bmx);
692: bmy=user->bmy, bmx=user->bmx;
694: info = DMDAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
695: info = VecSet(XL, lb); CHKERRQ(info);
696: info = DMDAVecGetArray(user->da,XL,&xl);CHKERRQ(info);
697: info = VecSet(XU, ub); CHKERRQ(info);
699: info = PetscOptionsHasName(PETSC_NULL,"-cylinder",&cylinder); CHKERRQ(info);
700: /* Compute the optional lower box */
701: if (cylinder){
702: for (i=xs; i< xs+xm; i++){
703: for (j=ys; j<ys+ym; j++){
704: t1=(2.0*i-mx)*bmy;
705: t2=(2.0*j-my)*bmx;
706: t3=bmx*bmx*bmy*bmy;
707: if ( t1*t1 + t2*t2 <= t3 ){
708: xl[j][i] = user->bheight;
709: }
710: }
711: }
712: } else {
713: /* Compute the optional lower box */
714: for (i=xs; i< xs+xm; i++){
715: for (j=ys; j<ys+ym; j++){
716: if (i>=(mx-bmx)/2 && i<mx-(mx-bmx)/2 &&
717: j>=(my-bmy)/2 && j<my-(my-bmy)/2 ){
718: xl[j][i] = user->bheight;
719: }
720: }
721: }
722: }
723:
724: info = DMDAVecRestoreArray(user->da,XL,&xl); CHKERRQ(info);
726: return(0);
727: }