Actual source code: ex16.c

  1: #include <petscsnes.h>
  2: #include <../src/snes/impls/vi/viimpl.h>
  3: #include <petscdmda.h>

  5: static  char help[]=
  6: "This example is an implementation of minimal surface area with \n\
  7: a plate problem from the TAO package (examples/plate2.c) \n\
  8: This example is based on a problem from the MINPACK-2 test suite.\n\
  9: Given a rectangular 2-D domain, boundary values along the edges of \n\
 10: the domain, and a plate represented by lower boundary conditions, \n\
 11: the objective is to find the surface with the minimal area that \n\
 12: satisfies the boundary conditions.\n\
 13: The command line options are:\n\
 14:   -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction\n\
 15:   -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction\n\
 16:   -bheight <ht>, where <ht> = height of the plate\n\
 17:   -start <st>, where <st> =0 for zero vector, <st> != 0 \n\
 18:                for an average of the boundary conditions\n\n";

 20: /*                                                                              
 21:    User-defined application context - contains data needed by the               
 22:    application-provided call-back routines, FormJacobian() and                  
 23:    FormFunction().                                                              
 24: */

 26: typedef struct {
 27:   DM           da;
 28:   Vec          Bottom, Top, Left, Right;
 29:   PetscScalar  bheight;
 30:   PetscInt     mx,my,bmx,bmy;
 31: } AppCtx;


 34: /* -------- User-defined Routines --------- */

 36: PetscErrorCode MSA_BoundaryConditions(AppCtx *);
 37: PetscErrorCode MSA_InitialPoint(AppCtx *, Vec);
 38: PetscErrorCode MSA_Plate(Vec,Vec,void*);
 39: PetscErrorCode FormGradient(SNES, Vec, Vec, void *);
 40: PetscErrorCode FormJacobian(SNES, Vec, Mat *, Mat*, MatStructure*,void *);

 44: int main(int argc, char **argv)
 45: {
 46:   PetscErrorCode  info;             /* used to check for functions returning nonzeros */
 47:   Vec             x,r;              /* solution and residual vectors */
 48:   Vec             xl,xu;            /* Bounds on the variables */
 49:   SNES            snes;             /* nonlinear solver context */
 50:   Mat             J;                /* Jacobian matrix */
 51:   PetscInt        N;            /* Number of elements in vector */
 52:   AppCtx          user;             /* user-defined work context */
 53:   PetscBool       flg;

 55:   /* Initialize PETSc */
 56:   PetscInitialize(&argc, &argv, (char *)0, help );

 58: #if defined(PETSC_USE_COMPLEX)
 59:   SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This example does not work for scalar type complex\n");
 60: #endif

 62:   /* Create distributed array to manage the 2d grid */
 63:   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);
 64:   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);CHKERRQ(info);

 66:   user.bheight=0.1;
 67:   info = PetscOptionsGetScalar(PETSC_NULL,"-bheight",&user.bheight,&flg); CHKERRQ(info);

 69:   user.bmx = user.mx/2; user.bmy = user.my/2;
 70:   info = PetscOptionsGetInt(PETSC_NULL,"-bmx",&user.bmx,&flg); CHKERRQ(info);
 71:   info = PetscOptionsGetInt(PETSC_NULL,"-bmy",&user.bmy,&flg); CHKERRQ(info);

 73:   PetscPrintf(PETSC_COMM_WORLD,"\n---- Minimum Surface Area With Plate Problem -----\n");
 74:   PetscPrintf(PETSC_COMM_WORLD,"mx:%d, my:%d, bmx:%d, bmy:%d, height:%4.2f\n",
 75:               user.mx,user.my,user.bmx,user.bmy,user.bheight);

 77:   /* Extract global vectors from DMDA; */
 78:   info = DMCreateGlobalVector(user.da,&x);CHKERRQ(info);
 79:   info = VecDuplicate(x, &r); CHKERRQ(info);

 81:   N = user.mx*user.my;
 82:   info = DMGetMatrix(user.da,MATAIJ,&J);CHKERRQ(info);

 84:   /* Create nonlinear solver context */
 85:   info = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(info);

 87:   /*  Set function evaluation and Jacobian evaluation  routines */
 88:   info = SNESSetFunction(snes,r,FormGradient,&user);CHKERRQ(info);
 89:   info = SNESSetJacobian(snes,J,J,FormJacobian,&user);CHKERRQ(info);

 91:   /* Set the boundary conditions */
 92:   info = MSA_BoundaryConditions(&user); CHKERRQ(info);

 94:   /* Set initial solution guess */
 95:   info = MSA_InitialPoint(&user, x); CHKERRQ(info);

 97:   info = SNESSetFromOptions(snes);CHKERRQ(info);

 99:   /* Set Bounds on variables */
100:   info = VecDuplicate(x, &xl); CHKERRQ(info);
101:   info = VecDuplicate(x, &xu); CHKERRQ(info);
102:   info = MSA_Plate(xl,xu,&user);CHKERRQ(info);

104:   info = SNESVISetVariableBounds(snes,xl,xu);CHKERRQ(info);

106:   /* Solve the application */
107:   info = SNESSolve(snes,PETSC_NULL,x);CHKERRQ(info);

109:   info = PetscOptionsHasName(PETSC_NULL,"-view_sol",&flg);CHKERRQ(info);
110:   if (flg) { info = VecView(x,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info); }

112:   /* Free memory */
113:   info = VecDestroy(&x); CHKERRQ(info);
114:   info = VecDestroy(&xl); CHKERRQ(info);
115:   info = VecDestroy(&xu); CHKERRQ(info);
116:   info = VecDestroy(&r); CHKERRQ(info);
117:   info = MatDestroy(&J); CHKERRQ(info);
118:   info = SNESDestroy(&snes); CHKERRQ(info);

120:   /* Free user-created data structures */
121:   info = DMDestroy(&user.da);CHKERRQ(info);
122:   info = VecDestroy(&user.Bottom); CHKERRQ(info);
123:   info = VecDestroy(&user.Top); CHKERRQ(info);
124:   info = VecDestroy(&user.Left); CHKERRQ(info);
125:   info = VecDestroy(&user.Right); CHKERRQ(info);

127:   info = PetscFinalize();

129:   return 0;
130: }

132: /* -------------------------------------------------------------------- */

136: /*  FormGradient - Evaluates gradient of f.             

138:     Input Parameters:
139: .   snes  - the SNES context
140: .   X     - input vector
141: .   ptr   - optional user-defined context, as set by SNESSetFunction()
142:     
143:     Output Parameters:
144: .   G - vector containing the newly evaluated gradient
145: */
146: PetscErrorCode FormGradient(SNES snes, Vec X, Vec G, void *ptr){
147:   AppCtx       *user = (AppCtx *) ptr;
148:   int          info;
149:   PetscInt     i,j;
150:   PetscInt     mx=user->mx, my=user->my;
151:   PetscScalar  hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
152:   PetscScalar  f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
153:   PetscScalar  df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
154:   PetscScalar  **g, **x;
155:   PetscInt     xs,xm,ys,ym;
156:   Vec          localX;
157:   PetscScalar  *top,*bottom,*left,*right;

160:   /* Initialize vector to zero */
161:   info = VecSet(G,0.0);CHKERRQ(info);

163:   /* Get local vector */
164:   info = DMGetLocalVector(user->da,&localX);CHKERRQ(info);
165:   info = VecGetArray(user->Top,&top); CHKERRQ(info);
166:   info = VecGetArray(user->Bottom,&bottom); CHKERRQ(info);
167:   info = VecGetArray(user->Left,&left); CHKERRQ(info);
168:   info = VecGetArray(user->Right,&right); CHKERRQ(info);

170:   /* Get ghost points */
171:   info = DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);CHKERRQ(info);
172:   info = DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);CHKERRQ(info);
173:   /* Get pointers to local vector data */
174:   info = DMDAVecGetArray(user->da,localX, &x); CHKERRQ(info);
175:   info = DMDAVecGetArray(user->da,G, &g); CHKERRQ(info);

177:   info = DMDAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);CHKERRQ(info);
178:   /* Compute function over the locally owned part of the mesh */
179:   for (j=ys; j < ys+ym; j++){
180:     for (i=xs; i< xs+xm; i++){
181: 
182:       xc = x[j][i];
183:       xlt=xrb=xl=xr=xb=xt=xc;
184: 
185:       if (i==0){ /* left side */
186:         xl= left[j-ys+1];
187:         xlt = left[j-ys+2];
188:       } else {
189:         xl = x[j][i-1];
190:       }

192:       if (j==0){ /* bottom side */
193:         xb=bottom[i-xs+1];
194:         xrb = bottom[i-xs+2];
195:       } else {
196:         xb = x[j-1][i];
197:       }
198: 
199:       if (i+1 == mx){ /* right side */
200:         xr=right[j-ys+1];
201:         xrb = right[j-ys];
202:       } else {
203:         xr = x[j][i+1];
204:       }

206:       if (j+1==0+my){ /* top side */
207:         xt=top[i-xs+1];
208:         xlt = top[i-xs];
209:       }else {
210:         xt = x[j+1][i];
211:       }

213:       if (i>0 && j+1<my){ /* left top side */
214:         xlt = x[j+1][i-1];
215:       }
216:       if (j>0 && i+1<mx){ /* right bottom */
217:         xrb = x[j-1][i+1];
218:       }

220:       d1 = (xc-xl);
221:       d2 = (xc-xr);
222:       d3 = (xc-xt);
223:       d4 = (xc-xb);
224:       d5 = (xr-xrb);
225:       d6 = (xrb-xb);
226:       d7 = (xlt-xl);
227:       d8 = (xt-xlt);
228: 
229:       df1dxc = d1*hydhx;
230:       df2dxc = ( d1*hydhx + d4*hxdhy );
231:       df3dxc = d3*hxdhy;
232:       df4dxc = ( d2*hydhx + d3*hxdhy );
233:       df5dxc = d2*hydhx;
234:       df6dxc = d4*hxdhy;

236:       d1 /= hx;
237:       d2 /= hx;
238:       d3 /= hy;
239:       d4 /= hy;
240:       d5 /= hy;
241:       d6 /= hx;
242:       d7 /= hy;
243:       d8 /= hx;

245:       f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
246:       f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
247:       f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
248:       f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
249:       f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
250:       f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);
251: 
252:       df1dxc /= f1;
253:       df2dxc /= f2;
254:       df3dxc /= f3;
255:       df4dxc /= f4;
256:       df5dxc /= f5;
257:       df6dxc /= f6;

259:       g[j][i] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc )/2.0;
260: 
261:     }
262:   }
263: 
264:   /* Restore vectors */
265:   info = DMDAVecRestoreArray(user->da,localX, &x); CHKERRQ(info);
266:   info = DMDAVecRestoreArray(user->da,G, &g); CHKERRQ(info);
267:   info = DMRestoreLocalVector(user->da,&localX);CHKERRQ(info);

269:   info = VecRestoreArray(user->Left,&left); CHKERRQ(info);
270:   info = VecRestoreArray(user->Top,&top); CHKERRQ(info);
271:   info = VecRestoreArray(user->Bottom,&bottom); CHKERRQ(info);
272:   info = VecRestoreArray(user->Right,&right); CHKERRQ(info);

274:   info = PetscLogFlops(67*mx*my); CHKERRQ(info);
275:   return(0);
276: }

278: /* ------------------------------------------------------------------- */
281: /*
282:    FormJacobian - Evaluates Jacobian matrix.

284:    Input Parameters:
285: .  snes - SNES context
286: .  X    - input vector
287: .  ptr  - optional user-defined context, as set by SNESSetJacobian()

289:    Output Parameters:
290: .  tH    - Jacobian matrix

292: */
293: PetscErrorCode FormJacobian(SNES snes, Vec X, Mat *tH, Mat* tHPre, MatStructure* flag, void *ptr)
294: {
295:   AppCtx          *user = (AppCtx *) ptr;
296:   Mat             H = *tH;
297:   PetscErrorCode  info;
298:   PetscInt        i,j,k;
299:   PetscInt        mx=user->mx, my=user->my;
300:   MatStencil      row,col[7];
301:   PetscScalar     hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
302:   PetscScalar     f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
303:   PetscScalar     hl,hr,ht,hb,hc,htl,hbr;
304:   PetscScalar     **x, v[7];
305:   PetscBool       assembled;
306:   PetscInt        xs,xm,ys,ym;
307:   Vec             localX;
308:   PetscScalar     *top,*bottom,*left,*right;

311:   /* Set various matrix options */
312:   info = MatAssembled(H,&assembled); CHKERRQ(info);
313:   if (assembled){info = MatZeroEntries(H);  CHKERRQ(info);}
314:   *flag=SAME_NONZERO_PATTERN;

316:   /* Get local vectors */
317:   info = DMGetLocalVector(user->da,&localX);CHKERRQ(info);
318:   info = VecGetArray(user->Top,&top); CHKERRQ(info);
319:   info = VecGetArray(user->Bottom,&bottom); CHKERRQ(info);
320:   info = VecGetArray(user->Left,&left); CHKERRQ(info);
321:   info = VecGetArray(user->Right,&right); CHKERRQ(info);

323:   /* Get ghost points */
324:   info = DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);CHKERRQ(info);
325:   info = DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);CHKERRQ(info);
326: 
327:   /* Get pointers to vector data */
328:   info = DMDAVecGetArray(user->da,localX, &x); CHKERRQ(info);

330:   info = DMDAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);CHKERRQ(info);
331:   /* Compute Jacobian over the locally owned part of the mesh */
332:   for (j=ys; j< ys+ym; j++){
333:     for (i=xs; i< xs+xm; i++){
334:       xc = x[j][i];
335:       xlt=xrb=xl=xr=xb=xt=xc;

337:       /* Left */
338:       if (i==0){
339:         xl= left[j+1];
340:         xlt = left[j+2];
341:       } else {
342:         xl = x[j][i-1];
343:       }
344: 
345:       /* Bottom */
346:       if (j==0){
347:         xb=bottom[i+1];
348:         xrb = bottom[i+2];
349:       } else {
350:         xb = x[j-1][i];
351:       }
352: 
353:       /* Right */
354:       if (i+1 == mx){
355:         xr=right[j+1];
356:         xrb = right[j];
357:       } else {
358:         xr = x[j][i+1];
359:       }

361:       /* Top */
362:       if (j+1==my){
363:         xt=top[i+1];
364:         xlt = top[i];
365:       }else {
366:         xt = x[j+1][i];
367:       }

369:       /* Top left */
370:       if (i>0 && j+1<my){
371:         xlt = x[j+1][i-1];
372:       }

374:       /* Bottom right */
375:       if (j>0 && i+1<mx){
376:         xrb = x[j-1][i+1];
377:       }

379:       d1 = (xc-xl)/hx;
380:       d2 = (xc-xr)/hx;
381:       d3 = (xc-xt)/hy;
382:       d4 = (xc-xb)/hy;
383:       d5 = (xrb-xr)/hy;
384:       d6 = (xrb-xb)/hx;
385:       d7 = (xlt-xl)/hy;
386:       d8 = (xlt-xt)/hx;
387: 
388:       f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
389:       f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
390:       f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
391:       f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
392:       f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
393:       f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);


396:       hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
397:         (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
398:       hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
399:         (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
400:       ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
401:         (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
402:       hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
403:         (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);

405:       hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
406:       htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);

408:       hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
409:         hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
410:         (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
411:         (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);

413:       hl/=2.0; hr/=2.0; ht/=2.0; hb/=2.0; hbr/=2.0; htl/=2.0;  hc/=2.0;

415:       k=0;
416:       row.i = i;row.j= j;
417:       /* Bottom */
418:       if (j>0){
419:         v[k]=hb;
420:         col[k].i = i; col[k].j=j-1; k++;
421:       }
422: 
423:       /* Bottom right */
424:       if (j>0 && i < mx -1){
425:         v[k]=hbr;
426:         col[k].i = i+1; col[k].j = j-1; k++;
427:       }
428: 
429:       /* left */
430:       if (i>0){
431:         v[k]= hl;
432:         col[k].i = i-1; col[k].j = j; k++;
433:       }
434: 
435:       /* Centre */
436:       v[k]= hc; col[k].i= row.i; col[k].j = row.j; k++;
437: 
438:       /* Right */
439:       if (i < mx-1 ){
440:         v[k]= hr;
441:         col[k].i= i+1; col[k].j = j;k++;
442:       }
443: 
444:       /* Top left */
445:       if (i>0 && j < my-1 ){
446:         v[k]= htl;
447:         col[k].i = i-1;col[k].j = j+1; k++;
448:       }
449: 
450:       /* Top */
451:       if (j < my-1 ){
452:         v[k]= ht;
453:         col[k].i = i; col[k].j = j+1; k++;
454:       }
455: 
456:       info = MatSetValuesStencil(H,1,&row,k,col,v,INSERT_VALUES);
457:       CHKERRQ(info);
458:     }
459:   }

461:   info = VecRestoreArray(user->Left,&left); CHKERRQ(info);
462:   info = VecRestoreArray(user->Top,&top); CHKERRQ(info);
463:   info = VecRestoreArray(user->Bottom,&bottom); CHKERRQ(info);
464:   info = VecRestoreArray(user->Right,&right); CHKERRQ(info);

466:   /* Assemble the matrix */
467:   info = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
468:   info = DMDAVecRestoreArray(user->da,localX,&x);CHKERRQ(info);
469:   info = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
470:   info = DMRestoreLocalVector(user->da,&localX);CHKERRQ(info);

472:   info = PetscLogFlops(199*mx*my); CHKERRQ(info);
473:   return(0);
474: }

476: /* ------------------------------------------------------------------- */
479: /* 
480:    MSA_BoundaryConditions -  Calculates the boundary conditions for
481:    the region.

483:    Input Parameter:
484: .  user - user-defined application context

486:    Output Parameter:
487: .  user - user-defined application context
488: */
489: PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
490: {
491:   PetscErrorCode  info;
492:   PetscInt        i,j,k,limit=0,maxits=5;
493:   PetscInt        mx=user->mx,my=user->my;
494:   PetscInt        xs,ys,xm,ym;
495:   PetscInt        bsize=0, lsize=0, tsize=0, rsize=0;
496:   PetscScalar     one=1.0, two=2.0, three=3.0, tol=1e-10;
497:   PetscScalar     fnorm,det,hx,hy,xt=0,yt=0;
498:   PetscScalar     u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
499:   PetscScalar     b=-0.5, t=0.5, l=-0.5, r=0.5;
500:   PetscScalar     *boundary;
501:   Vec             Bottom,Top,Right,Left;
502:   PetscScalar     scl=1.0;
503:   PetscBool       flg;


507:   info = DMDAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);

509:   bsize=xm+2; lsize=ym+2; rsize=ym+2; tsize=xm+2;

511:   info = VecCreateMPI(PETSC_COMM_WORLD,bsize,PETSC_DECIDE,&Bottom); CHKERRQ(info);
512:   info = VecCreateMPI(PETSC_COMM_WORLD,tsize,PETSC_DECIDE,&Top); CHKERRQ(info);
513:   info = VecCreateMPI(PETSC_COMM_WORLD,lsize,PETSC_DECIDE,&Left); CHKERRQ(info);
514:   info = VecCreateMPI(PETSC_COMM_WORLD,rsize,PETSC_DECIDE,&Right); CHKERRQ(info);

516:   user->Top=Top;
517:   user->Left=Left;
518:   user->Bottom=Bottom;
519:   user->Right=Right;

521:   hx= (r-l)/(mx+1); hy=(t-b)/(my+1);

523:   for (j=0; j<4; j++){
524:     if (j==0){
525:       yt=b;
526:       xt=l+hx*xs;
527:       limit=bsize;
528:       info = VecGetArray(Bottom,&boundary);CHKERRQ(info);
529:     } else if (j==1){
530:       yt=t;
531:       xt=l+hx*xs;
532:       limit=tsize;
533:       info = VecGetArray(Top,&boundary);CHKERRQ(info);
534:     } else if (j==2){
535:       yt=b+hy*ys;
536:       xt=l;
537:       limit=lsize;
538:       info = VecGetArray(Left,&boundary); CHKERRQ(info);
539:     } else { // if  (j==3)
540:       yt=b+hy*ys;
541:       xt=r;
542:       limit=rsize;
543:       info = VecGetArray(Right,&boundary);CHKERRQ(info);
544:     }

546:     for (i=0; i<limit; i++){
547:       u1=xt;
548:       u2=-yt;
549:       for (k=0; k<maxits; k++){
550:         nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
551:         nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
552:         fnorm=PetscSqrtReal(nf1*nf1+nf2*nf2);
553:         if (fnorm <= tol) break;
554:         njac11=one+u2*u2-u1*u1;
555:         njac12=two*u1*u2;
556:         njac21=-two*u1*u2;
557:         njac22=-one - u1*u1 + u2*u2;
558:         det = njac11*njac22-njac21*njac12;
559:         u1 = u1-(njac22*nf1-njac12*nf2)/det;
560:         u2 = u2-(njac11*nf2-njac21*nf1)/det;
561:       }

563:       boundary[i]=u1*u1-u2*u2;
564:       if (j==0 || j==1) {
565:         xt=xt+hx;
566:       } else { // if (j==2 || j==3)
567:         yt=yt+hy;
568:       }
569:     }

571:     if (j==0){
572:       info = VecRestoreArray(Bottom,&boundary); CHKERRQ(info);
573:     } else if (j==1){
574:       info = VecRestoreArray(Top,&boundary); CHKERRQ(info);
575:     } else if (j==2){
576:       info = VecRestoreArray(Left,&boundary); CHKERRQ(info);
577:     } else if (j==3){
578:       info = VecRestoreArray(Right,&boundary); CHKERRQ(info);
579:     }

581:   }

583:   /* Scale the boundary if desired */

585:   info = PetscOptionsGetReal(PETSC_NULL,"-bottom",&scl,&flg);
586:   CHKERRQ(info);
587:   if (flg){
588:     info = VecScale(Bottom, scl); CHKERRQ(info);
589:   }
590: 
591:   info = PetscOptionsGetReal(PETSC_NULL,"-top",&scl,&flg);
592:   CHKERRQ(info);
593:   if (flg){
594:     info = VecScale(Top, scl); CHKERRQ(info);
595:   }
596: 
597:   info = PetscOptionsGetReal(PETSC_NULL,"-right",&scl,&flg);
598:   CHKERRQ(info);
599:   if (flg){
600:     info = VecScale(Right, scl); CHKERRQ(info);
601:   }
602: 
603:   info = PetscOptionsGetReal(PETSC_NULL,"-left",&scl,&flg);
604:   CHKERRQ(info);
605:   if (flg){
606:     info = VecScale(Left, scl); CHKERRQ(info);
607:   }

609:   return(0);
610: }

612: /* ------------------------------------------------------------------- */
615: /*
616:    MSA_InitialPoint - Calculates the initial guess in one of three ways. 

618:    Input Parameters:
619: .  user - user-defined application context
620: .  X - vector for initial guess

622:    Output Parameters:
623: .  X - newly computed initial guess
624: */
625: PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
626: {
627:   PetscErrorCode  info;
628:   PetscInt        start=-1,i,j;
629:   PetscScalar     zero=0.0;
630:   PetscBool       flg;
631:   PetscScalar     *left,*right,*bottom,*top;

634:   info = PetscOptionsGetInt(PETSC_NULL,"-start",&start,&flg); CHKERRQ(info);

636:   if (flg && start==0){ /* The zero vector is reasonable */
637: 
638:     info = VecSet(X, zero); CHKERRQ(info);
639:     /* PLogInfo(user,"Min. Surface Area Problem: Start with 0 vector \n"); */


642:   } else { /* Take an average of the boundary conditions */
643:     PetscInt     mx=user->mx,my=user->my;
644:     PetscScalar  **x;
645:     PetscInt    xs,xm,ys,ym;
646: 
647:     info = VecGetArray(user->Top,&top); CHKERRQ(info);
648:     info = VecGetArray(user->Bottom,&bottom); CHKERRQ(info);
649:     info = VecGetArray(user->Left,&left); CHKERRQ(info);
650:     info = VecGetArray(user->Right,&right); CHKERRQ(info);

652:     /* Get pointers to vector data */
653:     info = DMDAVecGetArray(user->da,X,&x); CHKERRQ(info);
654:     info = DMDAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);CHKERRQ(info);

656:     /* Perform local computations */
657:     for (j=ys; j<ys+ym; j++){
658:       for (i=xs; i< xs+xm; i++){
659:         x[j][i] = ( (j+1)*bottom[i-xs+1]/my+(my-j+1)*top[i-xs+1]/(my+2)+
660:                    (i+1)*left[j-ys+1]/mx+(mx-i+1)*right[j-ys+1]/(mx+2))/2.0;
661:       }
662:     }
663: 
664:     /* Restore vectors */
665:     info = DMDAVecRestoreArray(user->da,X,&x); CHKERRQ(info);
666:     info = VecRestoreArray(user->Left,&left); CHKERRQ(info);
667:     info = VecRestoreArray(user->Top,&top); CHKERRQ(info);
668:     info = VecRestoreArray(user->Bottom,&bottom); CHKERRQ(info);
669:     info = VecRestoreArray(user->Right,&right); CHKERRQ(info);

671:   }
672:   return(0);
673: }

677: /* 
678:    MSA_Plate -  Calculates an obstacle for surface to stretch over.
679: */
680: PetscErrorCode MSA_Plate(Vec XL,Vec XU,void *ctx)
681: {
682:   AppCtx         *user=(AppCtx *)ctx;
683:   PetscErrorCode info;
684:   PetscInt       i,j;
685:   PetscInt       xs,ys,xm,ym;
686:   PetscInt       mx=user->mx, my=user->my, bmy, bmx;
687:   PetscScalar    t1,t2,t3;
688:   PetscScalar    **xl;
689:   PetscScalar    lb=-SNES_VI_INF, ub=SNES_VI_INF;
690:   PetscBool      cylinder;

692:   user->bmy = PetscMax(0,user->bmy);user->bmy = PetscMin(my,user->bmy);
693:   user->bmx = PetscMax(0,user->bmx);user->bmx = PetscMin(mx,user->bmx);
694:   bmy=user->bmy, bmx=user->bmx;

696:   info = DMDAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
697:   info = VecSet(XL, lb); CHKERRQ(info);
698:   info = DMDAVecGetArray(user->da,XL,&xl);CHKERRQ(info);
699:   info = VecSet(XU, ub); CHKERRQ(info);

701:   info = PetscOptionsHasName(PETSC_NULL,"-cylinder",&cylinder); CHKERRQ(info);
702:   /* Compute the optional lower box */
703:   if (cylinder){
704:     for (i=xs; i< xs+xm; i++){
705:       for (j=ys; j<ys+ym; j++){
706:         t1=(2.0*i-mx)*bmy;
707:         t2=(2.0*j-my)*bmx;
708:         t3=bmx*bmx*bmy*bmy;
709:         if ( t1*t1 + t2*t2 <= t3 ){
710:           xl[j][i] = user->bheight;
711:         }
712:       }
713:     }
714:   } else {
715:     /* Compute the optional lower box */
716:     for (i=xs; i< xs+xm; i++){
717:       for (j=ys; j<ys+ym; j++){
718:         if (i>=(mx-bmx)/2 && i<mx-(mx-bmx)/2 &&
719:             j>=(my-bmy)/2 && j<my-(my-bmy)/2 ){
720:           xl[j][i] = user->bheight;
721:         }
722:       }
723:     }
724:   }
725: 
726:   info = DMDAVecRestoreArray(user->da,XL,&xl); CHKERRQ(info);

728:   return(0);
729: }