Actual source code: plate2.c

petsc-3.8.3 2017-12-09
Report Typos and Errors
```  1:  #include <petscdmda.h>
2:  #include <petsctao.h>

4: static  char help[] =
5: "This example demonstrates use of the TAO package to \n\
6: solve an unconstrained minimization problem.  This example is based on a \n\
7: problem from the MINPACK-2 test suite.  Given a rectangular 2-D domain, \n\
8: boundary values along the edges of the domain, and a plate represented by \n\
9: lower boundary conditions, the objective is to find the\n\
10: surface with the minimal area that satisfies the boundary conditions.\n\
11: The command line options are:\n\
12:   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
13:   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\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 for random start, and <st> <0 \n\
18:                for an average of the boundary conditions\n\n";

20: /*T
21:    Concepts: TAO^Solving a bound constrained minimization problem
22:    Routines: TaoCreate();
24:    Routines: TaoSetHessianRoutine();
25:    Routines: TaoSetInitialVector();
26:    Routines: TaoSetVariableBounds();
27:    Routines: TaoSetFromOptions();
28:    Routines: TaoSolve(); TaoView();
29:    Routines: TaoDestroy();
30:    Processors: n
31: T*/

34: /*
35:    User-defined application context - contains data needed by the
37:    FormHessian().
38: */
39: typedef struct {
40:   /* problem parameters */
41:   PetscReal      bheight;                  /* Height of plate under the surface */
42:   PetscInt       mx, my;                   /* discretization in x, y directions */
43:   PetscInt       bmx,bmy;                  /* Size of plate under the surface */
44:   Vec            Bottom, Top, Left, Right; /* boundary values */

46:   /* Working space */
47:   Vec         localX, localV;           /* ghosted local vector */
48:   DM          dm;                       /* distributed array data structure */
49:   Mat         H;
50: } AppCtx;

52: /* -------- User-defined Routines --------- */

54: static PetscErrorCode MSA_BoundaryConditions(AppCtx*);
55: static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec);
56: static PetscErrorCode MSA_Plate(Vec,Vec,void*);
58: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);

60: /* For testing matrix free submatrices */
61: PetscErrorCode MatrixFreeHessian(Tao,Vec,Mat, Mat,void*);
62: PetscErrorCode MyMatMult(Mat,Vec,Vec);

64: int main( int argc, char **argv )
65: {
66:   PetscErrorCode         ierr;                 /* used to check for functions returning nonzeros */
67:   PetscInt               Nx, Ny;               /* number of processors in x- and y- directions */
68:   PetscInt               m, N;                 /* number of local and global elements in vectors */
69:   Vec                    x,xl,xu;               /* solution vector  and bounds*/
70:   PetscBool              flg;                /* A return variable when checking for user options */
71:   Tao                    tao;                  /* Tao solver context */
72:   ISLocalToGlobalMapping isltog;   /* local-to-global mapping object */
73:   Mat                    H_shell;                  /* to test matrix-free submatrices */
74:   AppCtx                 user;                 /* user-defined work context */

76:   /* Initialize PETSc, TAO */
77:   PetscInitialize( &argc, &argv,(char *)0,help );if (ierr) return ierr;

79:   /* Specify default dimension of the problem */
80:   user.mx = 10; user.my = 10; user.bheight=0.1;

82:   /* Check for any command line arguments that override defaults */
83:   PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,&flg);
84:   PetscOptionsGetInt(NULL,NULL,"-my",&user.my,&flg);
85:   PetscOptionsGetReal(NULL,NULL,"-bheight",&user.bheight,&flg);

87:   user.bmx = user.mx/2; user.bmy = user.my/2;
88:   PetscOptionsGetInt(NULL,NULL,"-bmx",&user.bmx,&flg);
89:   PetscOptionsGetInt(NULL,NULL,"-bmy",&user.bmy,&flg);

91:   PetscPrintf(PETSC_COMM_WORLD,"\n---- Minimum Surface Area With Plate Problem -----\n");
92:   PetscPrintf(PETSC_COMM_WORLD,"mx:%D, my:%D, bmx:%D, bmy:%D, height:%g\n",user.mx,user.my,user.bmx,user.bmy,(double)user.bheight);

94:   /* Calculate any derived values from parameters */
95:   N    = user.mx*user.my;

97:   /* Let Petsc determine the dimensions of the local vectors */
98:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

100:   /*
101:      A two dimensional distributed array will help define this problem,
102:      which derives from an elliptic PDE on two dimensional domain.  From
103:      the distributed array, Create the vectors.
104:   */
105:   DMDACreate2d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,user.mx,user.my,Nx,Ny,1,1,NULL,NULL,&user.dm);
106:   DMSetFromOptions(user.dm);
107:   DMSetUp(user.dm);
108:   /*
109:      Extract global and local vectors from DM; The local vectors are
110:      used solely as work space for the evaluation of the function,
111:      gradient, and Hessian.  Duplicate for remaining vectors that are
112:      the same types.
113:   */
114:   DMCreateGlobalVector(user.dm,&x); /* Solution */
115:   DMCreateLocalVector(user.dm,&user.localX);
116:   VecDuplicate(user.localX,&user.localV);

118:   VecDuplicate(x,&xl);
119:   VecDuplicate(x,&xu);

121:   /* The TAO code begins here */

123:   /*
124:      Create TAO solver and set desired solution method
125:      The method must either be TAOTRON or TAOBLMVM
126:      If TAOBLMVM is used, then hessian function is not called.
127:   */
128:   TaoCreate(PETSC_COMM_WORLD,&tao);
129:   TaoSetType(tao,TAOBLMVM);

131:   /* Set initial solution guess; */
132:   MSA_BoundaryConditions(&user);
133:   MSA_InitialPoint(&user,x);
134:   TaoSetInitialVector(tao,x);

136:   /* Set routines for function, gradient and hessian evaluation */

139:   VecGetLocalSize(x,&m);
140:   MatCreateAIJ(MPI_COMM_WORLD,m,m,N,N,7,NULL,3,NULL,&(user.H));
141:   MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);

143:   DMGetLocalToGlobalMapping(user.dm,&isltog);
144:   MatSetLocalToGlobalMapping(user.H,isltog,isltog);
145:   PetscOptionsHasName(NULL,NULL,"-matrixfree",&flg);
146:   if (flg) {
147:       MatCreateShell(PETSC_COMM_WORLD,m,m,N,N,(void*)&user,&H_shell);
148:       MatShellSetOperation(H_shell,MATOP_MULT,(void(*)(void))MyMatMult);
149:       MatSetOption(H_shell,MAT_SYMMETRIC,PETSC_TRUE);
150:       TaoSetHessianRoutine(tao,H_shell,H_shell,MatrixFreeHessian,(void*)&user);
151:   } else {
152:       TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);
153:   }

155:   /* Set Variable bounds */
156:   MSA_Plate(xl,xu,(void*)&user);
157:   TaoSetVariableBounds(tao,xl,xu);

159:   /* Check for any tao command line options */
160:   TaoSetFromOptions(tao);

162:   /* SOLVE THE APPLICATION */
163:   TaoSolve(tao);

165:   TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);

167:   /* Free TAO data structures */
168:   TaoDestroy(&tao);

170:   /* Free PETSc data structures */
171:   VecDestroy(&x);
172:   VecDestroy(&xl);
173:   VecDestroy(&xu);
174:   MatDestroy(&user.H);
175:   VecDestroy(&user.localX);
176:   VecDestroy(&user.localV);
177:   VecDestroy(&user.Bottom);
178:   VecDestroy(&user.Top);
179:   VecDestroy(&user.Left);
180:   VecDestroy(&user.Right);
181:   DMDestroy(&user.dm);
182:   if (flg) {
183:     MatDestroy(&H_shell);
184:   }
185:   PetscFinalize();
186:   return ierr;
187: }

191:     Input Parameters:
192: .   tao     - the Tao context
193: .   X      - input vector
194: .   userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()

196:     Output Parameters:
197: .   fcn     - the function value
198: .   G      - vector containing the newly evaluated gradient

200:    Notes:
201:    In this case, we discretize the domain and Create triangles. The
202:    surface of each triangle is planar, whose surface area can be easily
203:    computed.  The total surface area is found by sweeping through the grid
204:    and computing the surface area of the two triangles that have their
205:    right angle at the grid point.  The diagonal line segments on the
206:    grid that define the triangles run from top left to lower right.
207:    The numbering of points starts at the lower left and runs left to
208:    right, then bottom to top.
209: */
210: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G,void *userCtx)
211: {
212:   AppCtx         *user = (AppCtx *) userCtx;
214:   PetscInt       i,j,row;
215:   PetscInt       mx=user->mx, my=user->my;
216:   PetscInt       xs,xm,gxs,gxm,ys,ym,gys,gym;
217:   PetscReal      ft=0;
218:   PetscReal      zero=0.0;
219:   PetscReal      hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
220:   PetscReal      rhx=mx+1, rhy=my+1;
221:   PetscReal      f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
222:   PetscReal      df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
223:   PetscReal      *g, *x,*left,*right,*bottom,*top;
224:   Vec            localX = user->localX, localG = user->localV;

226:   /* Get local mesh boundaries */
227:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
228:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

230:   /* Scatter ghost points to local vector */
231:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
232:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

234:   /* Initialize vector to zero */
235:   VecSet(localG, zero);

237:   /* Get pointers to vector data */
238:   VecGetArray(localX,&x);
239:   VecGetArray(localG,&g);
240:   VecGetArray(user->Top,&top);
241:   VecGetArray(user->Bottom,&bottom);
242:   VecGetArray(user->Left,&left);
243:   VecGetArray(user->Right,&right);

245:   /* Compute function over the locally owned part of the mesh */
246:   for (j=ys; j<ys+ym; j++){
247:     for (i=xs; i< xs+xm; i++){
248:       row=(j-gys)*gxm + (i-gxs);

250:       xc = x[row];
251:       xlt=xrb=xl=xr=xb=xt=xc;

253:       if (i==0){ /* left side */
254:         xl= left[j-ys+1];
255:         xlt = left[j-ys+2];
256:       } else {
257:         xl = x[row-1];
258:       }

260:       if (j==0){ /* bottom side */
261:         xb=bottom[i-xs+1];
262:         xrb = bottom[i-xs+2];
263:       } else {
264:         xb = x[row-gxm];
265:       }

267:       if (i+1 == gxs+gxm){ /* right side */
268:         xr=right[j-ys+1];
269:         xrb = right[j-ys];
270:       } else {
271:         xr = x[row+1];
272:       }

274:       if (j+1==gys+gym){ /* top side */
275:         xt=top[i-xs+1];
276:         xlt = top[i-xs];
277:       }else {
278:         xt = x[row+gxm];
279:       }

281:       if (i>gxs && j+1<gys+gym){
282:         xlt = x[row-1+gxm];
283:       }
284:       if (j>gys && i+1<gxs+gxm){
285:         xrb = x[row+1-gxm];
286:       }

288:       d1 = (xc-xl);
289:       d2 = (xc-xr);
290:       d3 = (xc-xt);
291:       d4 = (xc-xb);
292:       d5 = (xr-xrb);
293:       d6 = (xrb-xb);
294:       d7 = (xlt-xl);
295:       d8 = (xt-xlt);

297:       df1dxc = d1*hydhx;
298:       df2dxc = ( d1*hydhx + d4*hxdhy );
299:       df3dxc = d3*hxdhy;
300:       df4dxc = ( d2*hydhx + d3*hxdhy );
301:       df5dxc = d2*hydhx;
302:       df6dxc = d4*hxdhy;

304:       d1 *= rhx;
305:       d2 *= rhx;
306:       d3 *= rhy;
307:       d4 *= rhy;
308:       d5 *= rhy;
309:       d6 *= rhx;
310:       d7 *= rhy;
311:       d8 *= rhx;

313:       f1 = PetscSqrtScalar( 1.0 + d1*d1 + d7*d7);
314:       f2 = PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
315:       f3 = PetscSqrtScalar( 1.0 + d3*d3 + d8*d8);
316:       f4 = PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
317:       f5 = PetscSqrtScalar( 1.0 + d2*d2 + d5*d5);
318:       f6 = PetscSqrtScalar( 1.0 + d4*d4 + d6*d6);

320:       ft = ft + (f2 + f4);

322:       df1dxc /= f1;
323:       df2dxc /= f2;
324:       df3dxc /= f3;
325:       df4dxc /= f4;
326:       df5dxc /= f5;
327:       df6dxc /= f6;

329:       g[row] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc ) * 0.5;

331:     }
332:   }

335:   /* Compute triangular areas along the border of the domain. */
336:   if (xs==0){ /* left side */
337:     for (j=ys; j<ys+ym; j++){
338:       d3=(left[j-ys+1] - left[j-ys+2])*rhy;
339:       d2=(left[j-ys+1] - x[(j-gys)*gxm])*rhx;
340:       ft = ft+PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
341:     }
342:   }
343:   if (ys==0){ /* bottom side */
344:     for (i=xs; i<xs+xm; i++){
345:       d2=(bottom[i+1-xs]-bottom[i-xs+2])*rhx;
346:       d3=(bottom[i-xs+1]-x[i-gxs])*rhy;
347:       ft = ft+PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
348:     }
349:   }

351:   if (xs+xm==mx){ /* right side */
352:     for (j=ys; j< ys+ym; j++){
353:       d1=(x[(j+1-gys)*gxm-1]-right[j-ys+1])*rhx;
354:       d4=(right[j-ys]-right[j-ys+1])*rhy;
355:       ft = ft+PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
356:     }
357:   }
358:   if (ys+ym==my){ /* top side */
359:     for (i=xs; i<xs+xm; i++){
360:       d1=(x[(gym-1)*gxm + i-gxs] - top[i-xs+1])*rhy;
361:       d4=(top[i-xs+1] - top[i-xs])*rhx;
362:       ft = ft+PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
363:     }
364:   }

366:   if (ys==0 && xs==0){
367:     d1=(left[0]-left[1])*rhy;
368:     d2=(bottom[0]-bottom[1])*rhx;
369:     ft +=PetscSqrtScalar( 1.0 + d1*d1 + d2*d2);
370:   }
371:   if (ys+ym == my && xs+xm == mx){
372:     d1=(right[ym+1] - right[ym])*rhy;
373:     d2=(top[xm+1] - top[xm])*rhx;
374:     ft +=PetscSqrtScalar( 1.0 + d1*d1 + d2*d2);
375:   }

377:   ft=ft*area;
378:   MPI_Allreduce(&ft,fcn,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);

381:   /* Restore vectors */
382:   VecRestoreArray(localX,&x);
383:   VecRestoreArray(localG,&g);
384:   VecRestoreArray(user->Left,&left);
385:   VecRestoreArray(user->Top,&top);
386:   VecRestoreArray(user->Bottom,&bottom);
387:   VecRestoreArray(user->Right,&right);

389:   /* Scatter values to global vector */
390:   DMLocalToGlobalBegin(user->dm,localG,INSERT_VALUES,G);
391:   DMLocalToGlobalEnd(user->dm,localG,INSERT_VALUES,G);

393:   PetscLogFlops(70*xm*ym);

395:   return 0;
396: }

398: /* ------------------------------------------------------------------- */
399: /*
400:    FormHessian - Evaluates Hessian matrix.

402:    Input Parameters:
403: .  tao  - the Tao context
404: .  x    - input vector
405: .  ptr  - optional user-defined context, as set by TaoSetHessianRoutine()

407:    Output Parameters:
408: .  A    - Hessian matrix
409: .  B    - optionally different preconditioning matrix

411:    Notes:
412:    Due to mesh point reordering with DMs, we must always work
413:    with the local mesh points, and then transform them to the new
414:    global numbering with the local-to-global mapping.  We cannot work
415:    directly with the global numbers for the original uniprocessor mesh!

417:    Two methods are available for imposing this transformation
418:    when setting matrix entries:
419:      (A) MatSetValuesLocal(), using the local ordering (including
420:          ghost points!)
421:          - Do the following two steps once, before calling TaoSolve()
422:            - Use DMGetISLocalToGlobalMapping() to extract the
423:              local-to-global map from the DM
424:            - Associate this map with the matrix by calling
425:              MatSetLocalToGlobalMapping()
426:          - Then set matrix entries using the local ordering
427:            by calling MatSetValuesLocal()
428:      (B) MatSetValues(), using the global ordering
429:          - Use DMGetGlobalIndices() to extract the local-to-global map
430:          - Then apply this map explicitly yourself
431:          - Set matrix entries using the global ordering by calling
432:            MatSetValues()
433:    Option (A) seems cleaner/easier in many cases, and is the procedure
434:    used in this example.
435: */
436: PetscErrorCode FormHessian(Tao tao,Vec X,Mat Hptr, Mat Hessian, void *ptr)
437: {
439:   AppCtx         *user = (AppCtx *) ptr;
440:   PetscInt       i,j,k,row;
441:   PetscInt       mx=user->mx, my=user->my;
442:   PetscInt       xs,xm,gxs,gxm,ys,ym,gys,gym,col[7];
443:   PetscReal      hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
444:   PetscReal      rhx=mx+1, rhy=my+1;
445:   PetscReal      f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
446:   PetscReal      hl,hr,ht,hb,hc,htl,hbr;
447:   PetscReal      *x,*left,*right,*bottom,*top;
448:   PetscReal      v[7];
449:   Vec            localX = user->localX;
450:   PetscBool      assembled;

453:   /* Set various matrix options */
454:   MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);

456:   /* Initialize matrix entries to zero */
457:   MatAssembled(Hessian,&assembled);
458:   if (assembled){MatZeroEntries(Hessian);}

460:   /* Get local mesh boundaries */
461:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
462:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

464:   /* Scatter ghost points to local vector */
465:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
466:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

468:   /* Get pointers to vector data */
469:   VecGetArray(localX,&x);
470:   VecGetArray(user->Top,&top);
471:   VecGetArray(user->Bottom,&bottom);
472:   VecGetArray(user->Left,&left);
473:   VecGetArray(user->Right,&right);

475:   /* Compute Hessian over the locally owned part of the mesh */

477:   for (i=xs; i< xs+xm; i++){

479:     for (j=ys; j<ys+ym; j++){

481:       row=(j-gys)*gxm + (i-gxs);

483:       xc = x[row];
484:       xlt=xrb=xl=xr=xb=xt=xc;

486:       /* Left side */
487:       if (i==gxs){
488:         xl= left[j-ys+1];
489:         xlt = left[j-ys+2];
490:       } else {
491:         xl = x[row-1];
492:       }

494:       if (j==gys){
495:         xb=bottom[i-xs+1];
496:         xrb = bottom[i-xs+2];
497:       } else {
498:         xb = x[row-gxm];
499:       }

501:       if (i+1 == gxs+gxm){
502:         xr=right[j-ys+1];
503:         xrb = right[j-ys];
504:       } else {
505:         xr = x[row+1];
506:       }

508:       if (j+1==gys+gym){
509:         xt=top[i-xs+1];
510:         xlt = top[i-xs];
511:       }else {
512:         xt = x[row+gxm];
513:       }

515:       if (i>gxs && j+1<gys+gym){
516:         xlt = x[row-1+gxm];
517:       }
518:       if (j>gys && i+1<gxs+gxm){
519:         xrb = x[row+1-gxm];
520:       }

523:       d1 = (xc-xl)*rhx;
524:       d2 = (xc-xr)*rhx;
525:       d3 = (xc-xt)*rhy;
526:       d4 = (xc-xb)*rhy;
527:       d5 = (xrb-xr)*rhy;
528:       d6 = (xrb-xb)*rhx;
529:       d7 = (xlt-xl)*rhy;
530:       d8 = (xlt-xt)*rhx;

532:       f1 = PetscSqrtScalar( 1.0 + d1*d1 + d7*d7);
533:       f2 = PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
534:       f3 = PetscSqrtScalar( 1.0 + d3*d3 + d8*d8);
535:       f4 = PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
536:       f5 = PetscSqrtScalar( 1.0 + d2*d2 + d5*d5);
537:       f6 = PetscSqrtScalar( 1.0 + d4*d4 + d6*d6);

540:       hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
541:         (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
542:       hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
543:         (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
544:       ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
545:         (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
546:       hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
547:         (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);

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

552:       hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
553:         hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
554:         (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
555:         (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);

557:       hl*=0.5; hr*=0.5; ht*=0.5; hb*=0.5; hbr*=0.5; htl*=0.5;  hc*=0.5;

559:       k=0;
560:       if (j>0){
561:         v[k]=hb; col[k]=row - gxm; k++;
562:       }

564:       if (j>0 && i < mx -1){
565:         v[k]=hbr; col[k]=row - gxm+1; k++;
566:       }

568:       if (i>0){
569:         v[k]= hl; col[k]=row - 1; k++;
570:       }

572:       v[k]= hc; col[k]=row; k++;

574:       if (i < mx-1 ){
575:         v[k]= hr; col[k]=row+1; k++;
576:       }

578:       if (i>0 && j < my-1 ){
579:         v[k]= htl; col[k] = row+gxm-1; k++;
580:       }

582:       if (j < my-1 ){
583:         v[k]= ht; col[k] = row+gxm; k++;
584:       }

586:       /*
587:          Set matrix values using local numbering, which was defined
588:          earlier, in the main routine.
589:       */
590:       MatSetValuesLocal(Hessian,1,&row,k,col,v,INSERT_VALUES);

592:     }
593:   }

595:   /* Restore vectors */
596:   VecRestoreArray(localX,&x);
597:   VecRestoreArray(user->Left,&left);
598:   VecRestoreArray(user->Top,&top);
599:   VecRestoreArray(user->Bottom,&bottom);
600:   VecRestoreArray(user->Right,&right);

602:   /* Assemble the matrix */
603:   MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
604:   MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);

606:   PetscLogFlops(199*xm*ym);
607:   return 0;
608: }

610: /* ------------------------------------------------------------------- */
611: /*
612:    MSA_BoundaryConditions -  Calculates the boundary conditions for
613:    the region.

615:    Input Parameter:
616: .  user - user-defined application context

618:    Output Parameter:
619: .  user - user-defined application context
620: */
621: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
622: {
623:   int        ierr;
624:   PetscInt   i,j,k,maxits=5,limit=0;
625:   PetscInt   xs,ys,xm,ym,gxs,gys,gxm,gym;
626:   PetscInt   mx=user->mx,my=user->my;
627:   PetscInt   bsize=0, lsize=0, tsize=0, rsize=0;
628:   PetscReal  one=1.0, two=2.0, three=3.0, scl=1.0, tol=1e-10;
629:   PetscReal  fnorm,det,hx,hy,xt=0,yt=0;
630:   PetscReal  u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
631:   PetscReal  b=-0.5, t=0.5, l=-0.5, r=0.5;
632:   PetscReal  *boundary;
633:   PetscBool  flg;
634:   Vec        Bottom,Top,Right,Left;

636:   /* Get local mesh boundaries */
637:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
638:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

641:   bsize=xm+2;
642:   lsize=ym+2;
643:   rsize=ym+2;
644:   tsize=xm+2;

646:   VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,&Bottom);
647:   VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,&Top);
648:   VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,&Left);
649:   VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,&Right);

651:   user->Top=Top;
652:   user->Left=Left;
653:   user->Bottom=Bottom;
654:   user->Right=Right;

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

658:   for (j=0; j<4; j++){
659:     if (j==0){
660:       yt=b;
661:       xt=l+hx*xs;
662:       limit=bsize;
663:       VecGetArray(Bottom,&boundary);
664:     } else if (j==1){
665:       yt=t;
666:       xt=l+hx*xs;
667:       limit=tsize;
668:       VecGetArray(Top,&boundary);
669:     } else if (j==2){
670:       yt=b+hy*ys;
671:       xt=l;
672:       limit=lsize;
673:       VecGetArray(Left,&boundary);
674:     } else if (j==3){
675:       yt=b+hy*ys;
676:       xt=r;
677:       limit=rsize;
678:       VecGetArray(Right,&boundary);
679:     }

681:     for (i=0; i<limit; i++){
682:       u1=xt;
683:       u2=-yt;
684:       for (k=0; k<maxits; k++){
685:         nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
686:         nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
687:         fnorm=PetscSqrtScalar(nf1*nf1+nf2*nf2);
688:         if (fnorm <= tol) break;
689:         njac11=one+u2*u2-u1*u1;
690:         njac12=two*u1*u2;
691:         njac21=-two*u1*u2;
692:         njac22=-one - u1*u1 + u2*u2;
693:         det = njac11*njac22-njac21*njac12;
694:         u1 = u1-(njac22*nf1-njac12*nf2)/det;
695:         u2 = u2-(njac11*nf2-njac21*nf1)/det;
696:       }

698:       boundary[i]=u1*u1-u2*u2;
699:       if (j==0 || j==1) {
700:         xt=xt+hx;
701:       } else if (j==2 || j==3){
702:         yt=yt+hy;
703:       }
704:     }
705:     if (j==0){
706:       VecRestoreArray(Bottom,&boundary);
707:     } else if (j==1){
708:       VecRestoreArray(Top,&boundary);
709:     } else if (j==2){
710:       VecRestoreArray(Left,&boundary);
711:     } else if (j==3){
712:       VecRestoreArray(Right,&boundary);
713:     }
714:   }

716:   /* Scale the boundary if desired */

718:   PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg);
719:   if (flg){
720:     VecScale(Bottom, scl);
721:   }
722:   PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg);
723:   if (flg){
724:     VecScale(Top, scl);
725:   }
726:   PetscOptionsGetReal(NULL,NULL,"-right",&scl,&flg);
727:   if (flg){
728:     VecScale(Right, scl);
729:   }

731:   PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg);
732:   if (flg){
733:     VecScale(Left, scl);
734:   }
735:   return 0;
736: }

739: /* ------------------------------------------------------------------- */
740: /*
741:    MSA_Plate -  Calculates an obstacle for surface to stretch over.

743:    Input Parameter:
744: .  user - user-defined application context

746:    Output Parameter:
747: .  user - user-defined application context
748: */
749: static PetscErrorCode MSA_Plate(Vec XL,Vec XU,void *ctx){

751:   AppCtx         *user=(AppCtx *)ctx;
753:   PetscInt       i,j,row;
754:   PetscInt       xs,ys,xm,ym;
755:   PetscInt       mx=user->mx, my=user->my, bmy, bmx;
756:   PetscReal      t1,t2,t3;
757:   PetscReal      *xl, lb=PETSC_NINFINITY, ub=PETSC_INFINITY;
758:   PetscBool      cylinder;

760:   user->bmy = PetscMax(0,user->bmy);user->bmy = PetscMin(my,user->bmy);
761:   user->bmx = PetscMax(0,user->bmx);user->bmx = PetscMin(mx,user->bmx);
762:   bmy=user->bmy; bmx=user->bmx;

764:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);

766:   VecSet(XL, lb);
767:   VecSet(XU, ub);

769:   VecGetArray(XL,&xl);

771:   PetscOptionsHasName(NULL,NULL,"-cylinder",&cylinder);
772:   /* Compute the optional lower box */
773:   if (cylinder){
774:     for (i=xs; i< xs+xm; i++){
775:       for (j=ys; j<ys+ym; j++){
776:         row=(j-ys)*xm + (i-xs);
777:         t1=(2.0*i-mx)*bmy;
778:         t2=(2.0*j-my)*bmx;
779:         t3=bmx*bmx*bmy*bmy;
780:         if ( t1*t1 + t2*t2 <= t3 ){
781:           xl[row] = user->bheight;
782:         }
783:       }
784:     }
785:   } else {
786:     /* Compute the optional lower box */
787:     for (i=xs; i< xs+xm; i++){
788:       for (j=ys; j<ys+ym; j++){
789:         row=(j-ys)*xm + (i-xs);
790:         if (i>=(mx-bmx)/2 && i<mx-(mx-bmx)/2 &&
791:             j>=(my-bmy)/2 && j<my-(my-bmy)/2 ){
792:           xl[row] = user->bheight;
793:         }
794:       }
795:     }
796:   }
797:     VecRestoreArray(XL,&xl);

799:   return 0;
800: }

803: /* ------------------------------------------------------------------- */
804: /*
805:    MSA_InitialPoint - Calculates the initial guess in one of three ways.

807:    Input Parameters:
808: .  user - user-defined application context
809: .  X - vector for initial guess

811:    Output Parameters:
812: .  X - newly computed initial guess
813: */
814: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
815: {
817:   PetscInt       start=-1,i,j;
818:   PetscReal      zero=0.0;
819:   PetscBool      flg;

821:   PetscOptionsGetInt(NULL,NULL,"-start",&start,&flg);
822:   if (flg && start==0){ /* The zero vector is reasonable */
823:     VecSet(X, zero);
824:   } else if (flg && start>0){ /* Try a random start between -0.5 and 0.5 */
825:     PetscRandom rctx;  PetscReal np5=-0.5;

827:     PetscRandomCreate(MPI_COMM_WORLD,&rctx);
828:     for (i=0; i<start; i++){
829:       VecSetRandom(X, rctx);
830:     }
831:     PetscRandomDestroy(&rctx);
832:     VecShift(X, np5);

834:   } else { /* Take an average of the boundary conditions */

836:     PetscInt row,xs,xm,gxs,gxm,ys,ym,gys,gym;
837:     PetscInt mx=user->mx,my=user->my;
838:     PetscReal *x,*left,*right,*bottom,*top;
839:     Vec    localX = user->localX;

841:     /* Get local mesh boundaries */
842:     DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
843:     DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

845:     /* Get pointers to vector data */
846:     VecGetArray(user->Top,&top);
847:     VecGetArray(user->Bottom,&bottom);
848:     VecGetArray(user->Left,&left);
849:     VecGetArray(user->Right,&right);

851:     VecGetArray(localX,&x);
852:     /* Perform local computations */
853:     for (j=ys; j<ys+ym; j++){
854:       for (i=xs; i< xs+xm; i++){
855:         row=(j-gys)*gxm + (i-gxs);
856:         x[row] = ( (j+1)*bottom[i-xs+1]/my + (my-j+1)*top[i-xs+1]/(my+2)+
857:                    (i+1)*left[j-ys+1]/mx + (mx-i+1)*right[j-ys+1]/(mx+2))/2.0;
858:       }
859:     }

861:     /* Restore vectors */
862:     VecRestoreArray(localX,&x);

864:     VecRestoreArray(user->Left,&left);
865:     VecRestoreArray(user->Top,&top);
866:     VecRestoreArray(user->Bottom,&bottom);
867:     VecRestoreArray(user->Right,&right);

869:     /* Scatter values into global vector */
870:     DMLocalToGlobalBegin(user->dm,localX,INSERT_VALUES,X);
871:     DMLocalToGlobalEnd(user->dm,localX,INSERT_VALUES,X);

873:   }
874:   return 0;
875: }

877: /* For testing matrix free submatrices */
878: PetscErrorCode MatrixFreeHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ptr)
879: {
881:   AppCtx         *user = (AppCtx*)ptr;
883:   FormHessian(tao,x,user->H,user->H,ptr);
884:   return(0);
885: }
886: PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y)
887: {
889:   void           *ptr;
890:   AppCtx         *user;
892:   MatShellGetContext(H_shell,&ptr);
893:   user = (AppCtx*)ptr;
894:   MatMult(user->H,X,Y);
895:   return(0);
896: }

```