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: }