Actual source code: ex6.c

  2: static char help[] = "Demonstrates using 3 DA's to manage a slightly non-trivial grid";

 4:  #include petscda.h

  6: struct _p_FA {
  7:   MPI_Comm   comm[3];
  8:   PetscInt   xl[3],yl[3],ml[3],nl[3];    /* corners and sizes of local vector in DA */
  9:   PetscInt   xg[3],yg[3],mg[3],ng[3];    /* corners and sizes of global vector in DA */
 10:   PetscInt   offl[3],offg[3];            /* offset in local and global vector of region 1, 2 and 3 portions */
 11:   Vec        g,l;
 12:   VecScatter vscat;
 13:   PetscInt   p1,p2,r1,r2,r1g,r2g,sw;
 14: };
 15: typedef struct _p_FA *FA;

 17: typedef struct {
 18:   PetscScalar X;
 19:   PetscScalar Y;
 20: } Field;

 22: PetscErrorCode FAGetLocalCorners(FA fa,PetscInt j,PetscInt *x,PetscInt *y,PetscInt *m,PetscInt *n)
 23: {
 25:   if (fa->comm[j]) {
 26:     *x = fa->xl[j];
 27:     *y = fa->yl[j];
 28:     *m = fa->ml[j];
 29:     *n = fa->nl[j];
 30:   } else {
 31:     *x = *y = *m = *n = 0;
 32:   }
 33:   return(0);
 34: }

 36: PetscErrorCode FAGetGlobalCorners(FA fa,PetscInt j,PetscInt *x,PetscInt *y,PetscInt *m,PetscInt *n)
 37: {
 39:   if (fa->comm[j]) {
 40:     *x = fa->xg[j];
 41:     *y = fa->yg[j];
 42:     *m = fa->mg[j];
 43:     *n = fa->ng[j];
 44:   } else {
 45:     *x = *y = *m = *n = 0;
 46:   }
 47:   return(0);
 48: }

 50: PetscErrorCode FAGetLocalArray(FA fa,Vec v,PetscInt j,Field ***f)
 51: {
 53:   PetscScalar    *va;
 54:   PetscInt       i;
 55:   Field          **a;

 58:   if (fa->comm[j]) {
 59:     VecGetArray(v,&va);
 60:     PetscMalloc(fa->nl[j]*sizeof(Field*),&a);
 61:     for (i=0; i<fa->nl[j]; i++) (a)[i] = (Field*) (va + 2*fa->offl[j] + i*2*fa->ml[j] - 2*fa->xl[j]);
 62:     *f = a - fa->yl[j];
 63:     VecRestoreArray(v,&va);
 64:   } else {
 65:     *f = 0;
 66:   }
 67:   return(0);
 68: }

 70: PetscErrorCode FARestoreLocalArray(FA fa,Vec v,PetscInt j,Field ***f)
 71: {
 73:   void           *dummy;

 76:   if (fa->comm[j]) {
 77:     dummy = *f + fa->yl[j];
 78:     PetscFree(dummy);
 79:   }
 80:   return(0);
 81: }

 83: PetscErrorCode FAGetGlobalArray(FA fa,Vec v,PetscInt j,Field ***f)
 84: {
 86:   PetscScalar    *va;
 87:   PetscInt       i;
 88:   Field          **a;

 91:   if (fa->comm[j]) {
 92:     VecGetArray(v,&va);
 93:     PetscMalloc(fa->ng[j]*sizeof(Field*),&a);
 94:     for (i=0; i<fa->ng[j]; i++) (a)[i] = (Field*) (va + 2*fa->offg[j] + i*2*fa->mg[j] - 2*fa->xg[j]);
 95:     *f = a - fa->yg[j];
 96:     VecRestoreArray(v,&va);
 97:   } else {
 98:     *f = 0;
 99:   }
100:   return(0);
101: }

103: PetscErrorCode FARestoreGlobalArray(FA fa,Vec v,PetscInt j,Field ***f)
104: {
106:   void           *dummy;

109:   if (fa->comm[j]) {
110:     dummy = *f + fa->yg[j];
111:     PetscFree(dummy);
112:   }
113:   return(0);
114: }

116: PetscErrorCode FAGetGlobalVector(FA fa,Vec *v)
117: {
120:   VecDuplicate(fa->g,v);
121:   return(0);
122: }

124: PetscErrorCode FAGetLocalVector(FA fa,Vec *v)
125: {
128:   VecDuplicate(fa->l,v);
129:   return(0);
130: }

132: PetscErrorCode FAGlobalToLocal(FA fa,Vec g,Vec l)
133: {
136:   VecScatterBegin(fa->vscat,g,l,INSERT_VALUES,SCATTER_FORWARD);
137:   VecScatterEnd(fa->vscat,g,l,INSERT_VALUES,SCATTER_FORWARD);
138:   return(0);
139: }

141: PetscErrorCode FADestroy(FA fa)
142: {
145:   VecDestroy(fa->g);
146:   VecDestroy(fa->l);
147:   VecScatterDestroy(fa->vscat);
148:   PetscFree(fa);
149:   return(0);
150: }

152: PetscErrorCode FACreate(FA *infa)
153: {
154:   FA             fa;
155:   PetscMPIInt    rank;
156:   PetscInt       tonglobal,globalrstart,x,nx,y,ny,*tonatural,i,j,*to,*from,offt[3];
157:   PetscInt       *fromnatural,fromnglobal,nscat,nlocal,cntl1,cntl2,cntl3,*indices;

160:   /* Each DA manages the local vector for the portion of region 1, 2, and 3 for that processor
161:      Each DA can belong on any subset (overlapping between DA's or not) of processors
162:      For processes that a particular DA does not exist on, the corresponding comm should be set to zero
163:   */
164:   DA             da1 = 0,da2 = 0,da3 = 0;
165:   /* 
166:       v1, v2, v3 represent the local vector for a single DA
167:   */
168:   Vec            vl1 = 0,vl2 = 0,vl3 = 0, vg1 = 0, vg2 = 0,vg3 = 0;

170:   /*
171:      globalvec and friends represent the global vectors that are used for the PETSc solvers 
172:      localvec represents the concatenation of the (up to) 3 local vectors; vl1, vl2, vl3

174:      tovec and friends represent intermediate vectors that are ONLY used for setting up the 
175:      final communication patterns. Once this setup routine is complete they are destroyed.
176:      The tovec  is like the globalvec EXCEPT it has redundant locations for the ghost points
177:      between regions 2+3 and 1.
178:   */
179:   AO             toao,globalao;
180:   IS             tois,globalis,is;
181:   Vec            tovec,globalvec,localvec;
182:   VecScatter     vscat;
183:   PetscScalar    *globalarray,*localarray,*toarray;

185:   PetscNew(struct _p_FA,&fa);
186:   /*
187:       fa->sw is the stencil width  

189:       fa->p1 is the width of region 1, fa->p2 the width of region 2 (must be the same)
190:       fa->r1 height of region 1 
191:       fa->r2 height of region 2
192:  
193:       fa->r2 is also the height of region 3-4
194:       (fa->p1 - fa->p2)/2 is the width of both region 3 and region 4
195:   */
196:   fa->p1 = 24;
197:   fa->p2 = 15;
198:   fa->r1 = 6;
199:   fa->r2 = 6;
200:   fa->sw = 1;
201:   fa->r1g = fa->r1 + fa->sw;
202:   fa->r2g = fa->r2 + fa->sw;

204:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

206:   fa->comm[0] = PETSC_COMM_WORLD;
207:   fa->comm[1] = PETSC_COMM_WORLD;
208:   fa->comm[2] = PETSC_COMM_WORLD;
209:   /* Test case with different communicators */
210:   /* Normally one would use MPI_Comm routines to build MPI communicators on which you wish to partition the DAs*/
211:   /*
212:   if (rank == 0) {
213:     fa->comm[0] = PETSC_COMM_SELF;
214:     fa->comm[1] = 0;
215:     fa->comm[2] = 0;
216:   } else if (rank == 1) {
217:     fa->comm[0] = 0;  
218:     fa->comm[1] = PETSC_COMM_SELF;
219:     fa->comm[2] = 0;  
220:   } else {
221:     fa->comm[0] = 0;  
222:     fa->comm[1] = 0;
223:     fa->comm[2] = PETSC_COMM_SELF;
224:   } */

226:   if (fa->p2 > fa->p1 - 3)   SETERRQ(1,"Width of region fa->p2 must be at least 3 less then width of region 1");
227:   if (!((fa->p2 - fa->p1) % 2)) SETERRQ(1,"width of region 3 must NOT be divisible by 2");

229:   if (fa->comm[1]) {
230:     DACreate2d(fa->comm[1],DA_XPERIODIC,DA_STENCIL_BOX,fa->p2,fa->r2g,PETSC_DECIDE,PETSC_DECIDE,1,fa->sw,PETSC_NULL,PETSC_NULL,&da2);
231:     DAGetLocalVector(da2,&vl2);
232:     DAGetGlobalVector(da2,&vg2);
233:   }
234:   if (fa->comm[2]) {
235:     DACreate2d(fa->comm[2],DA_NONPERIODIC,DA_STENCIL_BOX,fa->p1-fa->p2,fa->r2g,PETSC_DECIDE,PETSC_DECIDE,1,fa->sw,PETSC_NULL,PETSC_NULL,&da3);
236:     DAGetLocalVector(da3,&vl3);
237:     DAGetGlobalVector(da3,&vg3);
238:   }
239:   if (fa->comm[0]) {
240:     DACreate2d(fa->comm[0],DA_NONPERIODIC,DA_STENCIL_BOX,fa->p1,fa->r1g,PETSC_DECIDE,PETSC_DECIDE,1,fa->sw,PETSC_NULL,PETSC_NULL,&da1);
241:     DAGetLocalVector(da1,&vl1);
242:     DAGetGlobalVector(da1,&vg1);
243:   }

245:   /* count the number of unknowns owned on each processor and determine the starting point of each processors ownership 
246:      for global vector with redundancy */
247:   tonglobal = 0;
248:   if (fa->comm[1]) {
249:     DAGetCorners(da2,&x,&y,0,&nx,&ny,0);
250:     tonglobal += nx*ny;
251:   }
252:   if (fa->comm[2]) {
253:     DAGetCorners(da3,&x,&y,0,&nx,&ny,0);
254:     tonglobal += nx*ny;
255:   }
256:   if (fa->comm[0]) {
257:     DAGetCorners(da1,&x,&y,0,&nx,&ny,0);
258:     tonglobal += nx*ny;
259:   }
260:   PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] Number of unknowns owned %d\n",rank,tonglobal);
261:   PetscSynchronizedFlush(PETSC_COMM_WORLD);
262: 
263:   /* Get tonatural number for each node */
264:   PetscMalloc((tonglobal+1)*sizeof(PetscInt),&tonatural);
265:   tonglobal = 0;
266:   if (fa->comm[1]) {
267:     DAGetCorners(da2,&x,&y,0,&nx,&ny,0);
268:     for (j=0; j<ny; j++) {
269:       for (i=0; i<nx; i++) {
270:         tonatural[tonglobal++] = (fa->p1 - fa->p2)/2 + x + i + fa->p1*(y + j);
271:       }
272:     }
273:   }
274:   if (fa->comm[2]) {
275:     DAGetCorners(da3,&x,&y,0,&nx,&ny,0);
276:     for (j=0; j<ny; j++) {
277:       for (i=0; i<nx; i++) {
278:         if (x + i < (fa->p1 - fa->p2)/2) tonatural[tonglobal++] = x + i + fa->p1*(y + j);
279:         else tonatural[tonglobal++] = fa->p2 + x + i + fa->p1*(y + j);
280:       }
281:     }
282:   }
283:   if (fa->comm[0]) {
284:     DAGetCorners(da1,&x,&y,0,&nx,&ny,0);
285:     for (j=0; j<ny; j++) {
286:       for (i=0; i<nx; i++) {
287:         tonatural[tonglobal++] = fa->p1*fa->r2g + x + i + fa->p1*(y + j);
288:       }
289:     }
290:   }
291:   /*  PetscIntView(tonglobal,tonatural,PETSC_VIEWER_STDOUT_WORLD); */
292:   AOCreateBasic(PETSC_COMM_WORLD,tonglobal,tonatural,0,&toao);
293:   PetscFree(tonatural);

295:   /* count the number of unknowns owned on each processor and determine the starting point of each processors ownership 
296:      for global vector without redundancy */
297:   fromnglobal = 0;
298:   fa->offg[1] = 0;
299:   offt[1]     = 0;
300:   if (fa->comm[1]) {
301:     DAGetCorners(da2,&x,&y,0,&nx,&ny,0);
302:     offt[2] = nx*ny;
303:     if (y+ny == fa->r2g) {ny--;}  /* includes the ghost points on the upper side */
304:     fromnglobal += nx*ny;
305:     fa->offg[2] = fromnglobal;
306:   } else {
307:     offt[2] = 0;
308:     fa->offg[2] = 0;
309:   }
310:   if (fa->comm[2]) {
311:     DAGetCorners(da3,&x,&y,0,&nx,&ny,0);
312:     offt[0] = offt[2] + nx*ny;
313:     if (y+ny == fa->r2g) {ny--;}  /* includes the ghost points on the upper side */
314:     fromnglobal += nx*ny;
315:     fa->offg[0] = fromnglobal;
316:   } else {
317:     offt[0]     = offt[2];
318:     fa->offg[0] = fromnglobal;
319:   }
320:   if (fa->comm[0]) {
321:     DAGetCorners(da1,&x,&y,0,&nx,&ny,0);
322:     if (y == 0) {ny--;}  /* includes the ghost points on the lower side */
323:     fromnglobal += nx*ny;
324:   }
325:   MPI_Scan(&fromnglobal,&globalrstart,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);
326:   globalrstart -= fromnglobal;
327:   PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] Number of unknowns owned %d\n",rank,fromnglobal);
328:   PetscSynchronizedFlush(PETSC_COMM_WORLD);

330:   /* Get fromnatural number for each node */
331:   PetscMalloc((fromnglobal+1)*sizeof(PetscInt),&fromnatural);
332:   fromnglobal = 0;
333:   if (fa->comm[1]) {
334:     DAGetCorners(da2,&x,&y,0,&nx,&ny,0);
335:     if (y+ny == fa->r2g) {ny--;}  /* includes the ghost points on the upper side */
336:     fa->xg[1] = x; fa->yg[1] = y; fa->mg[1] = nx; fa->ng[1] = ny;
337:     DAGetGhostCorners(da2,&fa->xl[1],&fa->yl[1],0,&fa->ml[1],&fa->nl[1],0);
338:     for (j=0; j<ny; j++) {
339:       for (i=0; i<nx; i++) {
340:         fromnatural[fromnglobal++] = (fa->p1 - fa->p2)/2 + x + i + fa->p1*(y + j);
341:       }
342:     }
343:   }
344:   if (fa->comm[2]) {
345:     DAGetCorners(da3,&x,&y,0,&nx,&ny,0);
346:     if (y+ny == fa->r2g) {ny--;}  /* includes the ghost points on the upper side */
347:     fa->xg[2] = x; fa->yg[2] = y; fa->mg[2] = nx; fa->ng[2] = ny;
348:     DAGetGhostCorners(da3,&fa->xl[2],&fa->yl[2],0,&fa->ml[2],&fa->nl[2],0);
349:     for (j=0; j<ny; j++) {
350:       for (i=0; i<nx; i++) {
351:         if (x + i < (fa->p1 - fa->p2)/2) fromnatural[fromnglobal++] = x + i + fa->p1*(y + j);
352:         else fromnatural[fromnglobal++] = fa->p2 + x + i + fa->p1*(y + j);
353:       }
354:     }
355:   }
356:   if (fa->comm[0]) {
357:     DAGetCorners(da1,&x,&y,0,&nx,&ny,0);
358:     if (y == 0) {ny--;}  /* includes the ghost points on the lower side */
359:     else y--;
360:     fa->xg[0] = x; fa->yg[0] = y; fa->mg[0] = nx; fa->ng[0] = ny;
361:     DAGetGhostCorners(da1,&fa->xl[0],&fa->yl[0],0,&fa->ml[0],&fa->nl[0],0);
362:     for (j=0; j<ny; j++) {
363:       for (i=0; i<nx; i++) {
364:         fromnatural[fromnglobal++] = fa->p1*fa->r2 + x + i + fa->p1*(y + j);
365:       }
366:     }
367:   }
368:   /*PetscIntView(fromnglobal,fromnatural,PETSC_VIEWER_STDOUT_WORLD);*/
369:   AOCreateBasic(PETSC_COMM_WORLD,fromnglobal,fromnatural,0,&globalao);
370:   PetscFree(fromnatural);

372:   /* ---------------------------------------------------*/
373:   /* Create the scatter that updates 1 from 2 and 3 and 3 and 2 from 1 */
374:   /* currently handles stencil width of 1 ONLY */
375:   PetscMalloc(tonglobal*sizeof(PetscInt),&to);
376:   PetscMalloc(tonglobal*sizeof(PetscInt),&from);
377:   nscat = 0;
378:   if (fa->comm[1]) {
379:     DAGetCorners(da2,&x,&y,0,&nx,&ny,0);
380:     for (j=0; j<ny; j++) {
381:       for (i=0; i<nx; i++) {
382:         to[nscat] = from[nscat] = (fa->p1 - fa->p2)/2 + x + i + fa->p1*(y + j);nscat++;
383:       }
384:     }
385:   }
386:   if (fa->comm[2]) {
387:     DAGetCorners(da3,&x,&y,0,&nx,&ny,0);
388:     for (j=0; j<ny; j++) {
389:       for (i=0; i<nx; i++) {
390:         if (x + i < (fa->p1 - fa->p2)/2) {
391:           to[nscat]   = from[nscat] = x + i + fa->p1*(y + j);nscat++;
392:         } else {
393:           to[nscat]   = from[nscat] = fa->p2 + x + i + fa->p1*(y + j);nscat++;
394:         }
395:       }
396:     }
397:   }
398:   if (fa->comm[0]) {
399:     DAGetCorners(da1,&x,&y,0,&nx,&ny,0);
400:     for (j=0; j<ny; j++) {
401:       for (i=0; i<nx; i++) {
402:         to[nscat]     = fa->p1*fa->r2g + x + i + fa->p1*(y + j);
403:         from[nscat++] = fa->p1*(fa->r2 - 1) + x + i + fa->p1*(y + j);
404:       }
405:     }
406:   }
407:   AOApplicationToPetsc(toao,nscat,to);
408:   AOApplicationToPetsc(globalao,nscat,from);
409:   ISCreateGeneral(PETSC_COMM_WORLD,nscat,to,&tois);
410:   ISCreateGeneral(PETSC_COMM_WORLD,nscat,from,&globalis);
411:   PetscFree(to);
412:   PetscFree(from);
413:   VecCreateMPI(PETSC_COMM_WORLD,tonglobal,PETSC_DETERMINE,&tovec);
414:   VecCreateMPI(PETSC_COMM_WORLD,fromnglobal,PETSC_DETERMINE,&globalvec);
415:   VecScatterCreate(globalvec,globalis,tovec,tois,&vscat);
416:   ISDestroy(tois);
417:   ISDestroy(globalis);
418:   AODestroy(globalao);
419:   AODestroy(toao);

421:   /* fill up global vector without redundant values with PETSc global numbering */
422:   VecGetArray(globalvec,&globalarray);
423:   for (i=0; i<fromnglobal; i++) {
424:     globalarray[i] = globalrstart + i;
425:   }
426:   VecRestoreArray(globalvec,&globalarray);
427: 
428:   /* scatter PETSc global indices to redundant valueed array */
429:   VecScatterBegin(vscat,globalvec,tovec,INSERT_VALUES,SCATTER_FORWARD);
430:   VecScatterEnd(vscat,globalvec,tovec,INSERT_VALUES,SCATTER_FORWARD);
431: 
432:   /* Create local vector that is the concatenation of the local vectors */
433:   nlocal = 0;
434:   cntl1  = cntl2 = cntl3 = 0;
435:   if (fa->comm[1]) {
436:     VecGetSize(vl2,&cntl2);
437:     nlocal += cntl2;
438:   }
439:   if (fa->comm[2]) {
440:     VecGetSize(vl3,&cntl3);
441:     nlocal += cntl3;
442:   }
443:   if (fa->comm[0]) {
444:     VecGetSize(vl1,&cntl1);
445:     nlocal += cntl1;
446:   }
447:   fa->offl[0] = cntl2 + cntl3;
448:   fa->offl[1] = 0;
449:   fa->offl[2] = cntl2;
450:   VecCreateSeq(PETSC_COMM_SELF,nlocal,&localvec);
451: 
452:   /* cheat so that  vl1, vl2, vl3 shared array memory with localvec */
453:   VecGetArray(localvec,&localarray);
454:   VecGetArray(tovec,&toarray);
455:   if (fa->comm[1]) {
456:     VecPlaceArray(vl2,localarray+fa->offl[1]);
457:     VecPlaceArray(vg2,toarray+offt[1]);
458:     DAGlobalToLocalBegin(da2,vg2,INSERT_VALUES,vl2);
459:     DAGlobalToLocalEnd(da2,vg2,INSERT_VALUES,vl2);
460:     DARestoreGlobalVector(da2,&vg2);
461:   }
462:   if (fa->comm[2]) {
463:     VecPlaceArray(vl3,localarray+fa->offl[2]);
464:     VecPlaceArray(vg3,toarray+offt[2]);
465:     DAGlobalToLocalBegin(da3,vg3,INSERT_VALUES,vl3);
466:     DAGlobalToLocalEnd(da3,vg3,INSERT_VALUES,vl3);
467:     DARestoreGlobalVector(da3,&vg3);
468:   }
469:   if (fa->comm[0]) {
470:     VecPlaceArray(vl1,localarray+fa->offl[0]);
471:     VecPlaceArray(vg1,toarray+offt[0]);
472:     DAGlobalToLocalBegin(da1,vg1,INSERT_VALUES,vl1);
473:     DAGlobalToLocalEnd(da1,vg1,INSERT_VALUES,vl1);
474:     DARestoreGlobalVector(da1,&vg1);
475:   }
476:   VecRestoreArray(localvec,&localarray);
477:   VecRestoreArray(tovec,&toarray);

479:   /* no longer need the redundant vector and VecScatter to it */
480:   VecScatterDestroy(vscat);
481:   VecDestroy(tovec);

483:   /* Create final scatter that goes directly from globalvec to localvec */
484:   /* this is the one to be used in the application code */
485:   PetscMalloc(nlocal*sizeof(PetscInt),&indices);
486:   VecGetArray(localvec,&localarray);
487:   for (i=0; i<nlocal; i++) {
488:     indices[i] = (PetscInt) (2*localarray[i]);
489:   }
490:   VecRestoreArray(localvec,&localarray);
491:   ISCreateBlock(PETSC_COMM_WORLD,2,nlocal,indices,&is);
492:   PetscFree(indices);

494:   VecCreateSeq(PETSC_COMM_SELF,2*nlocal,&fa->l);
495:   VecCreateMPI(PETSC_COMM_WORLD,2*fromnglobal,PETSC_DETERMINE,&fa->g);

497:   VecScatterCreate(fa->g,is,fa->l,PETSC_NULL,&fa->vscat);
498:   ISDestroy(is);

500:   VecDestroy(globalvec);
501:   VecDestroy(localvec);
502:   if (fa->comm[0]) {
503:     DARestoreLocalVector(da1,&vl1);
504:     DADestroy(da1);
505:   }
506:   if (fa->comm[1]) {
507:     DARestoreLocalVector(da2,&vl2);
508:     DADestroy(da2);
509:   }
510:   if (fa->comm[2]) {
511:     DARestoreLocalVector(da3,&vl3);
512:     DADestroy(da3);
513:   }
514:   *infa = fa;
515:   return(0);
516: }

518: /* Crude graphics to test that the ghost points are properly updated */
519:  #include petscdraw.h

521: typedef struct {
522:   PetscInt     m[3],n[3];
523:   PetscScalar  *xy[3];
524: } ZoomCtx;

526: PetscErrorCode DrawPatch(PetscDraw draw,void *ctx)
527: {
528:   ZoomCtx        *zctx = (ZoomCtx*)ctx;
530:   PetscInt       m,n,i,j,id,k;
531:   PetscReal      x1,x2,x3,x4,y_1,y2,y3,y4;
532:   PetscScalar    *xy;

535:   for (k=0; k<3; k++) {
536:     m    = zctx->m[k];
537:     n    = zctx->n[k];
538:     xy   = zctx->xy[k];

540:     for (j=0; j<n-1; j++) {
541:       for (i=0; i<m-1; i++) {
542:         id = i+j*m;    x1 = xy[2*id];y_1 = xy[2*id+1];
543:         id = i+j*m+1;  x2 = xy[2*id];y2  = xy[2*id+1];
544:         id = i+j*m+1+m;x3 = xy[2*id];y3  = xy[2*id+1];
545:         id = i+j*m+m;  x4 = xy[2*id];y4  = xy[2*id+1];
546:         PetscDrawLine(draw,x1,y_1,x2,y2,PETSC_DRAW_BLACK);
547:         PetscDrawLine(draw,x2,y2,x3,y3,PETSC_DRAW_BLACK);
548:         PetscDrawLine(draw,x3,y3,x4,y4,PETSC_DRAW_BLACK);
549:         PetscDrawLine(draw,x4,y4,x1,y_1,PETSC_DRAW_BLACK);
550:       }
551:     }
552:   }
553:   PetscDrawFlush(draw);
554:   return(0);
555: }

557: PetscErrorCode DrawFA(FA fa,Vec v)
558: {
560:   PetscScalar    *va;
561:   ZoomCtx        zctx;
562:   PetscDraw      draw;
563:   PetscReal      xmint = 10000.0,xmaxt = -10000.0,ymint = 100000.0,ymaxt = -10000.0;
564:   PetscReal      xmin,xmax,ymin,ymax;
565:   PetscInt       i,vn,ln,j;

568:   VecGetArray(v,&va);
569:   VecGetSize(v,&vn);
570:   VecGetSize(fa->l,&ln);
571:   for (j=0; j<3; j++) {
572:     if (vn == ln) {
573:       zctx.xy[j] = va + 2*fa->offl[j];
574:       zctx.m[j]  = fa->ml[j];
575:       zctx.n[j]  = fa->nl[j];
576:     } else {
577:       zctx.xy[j] = va + 2*fa->offg[j];
578:       zctx.m[j]  = fa->mg[j];
579:       zctx.n[j]  = fa->ng[j];
580:     }
581:     for (i=0; i<zctx.m[j]*zctx.n[j]; i++) {
582:       if (zctx.xy[j][2*i] > xmax) xmax = zctx.xy[j][2*i];
583:       if (zctx.xy[j][2*i] < xmin) xmin = zctx.xy[j][2*i];
584:       if (zctx.xy[j][2*i+1] > ymax) ymax = zctx.xy[j][2*i+1];
585:       if (zctx.xy[j][2*i+1] < ymin) ymin = zctx.xy[j][2*i+1];
586:     }
587:   }
588:   MPI_Allreduce(&xmin,&xmint,1,MPI_DOUBLE,MPI_MIN,PETSC_COMM_WORLD);
589:   MPI_Allreduce(&xmax,&xmaxt,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD);
590:   MPI_Allreduce(&ymin,&ymint,1,MPI_DOUBLE,MPI_MIN,PETSC_COMM_WORLD);
591:   MPI_Allreduce(&ymax,&ymaxt,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD);
592:   xmin = xmint - .2*(xmaxt - xmint);
593:   xmax = xmaxt + .2*(xmaxt - xmint);
594:   ymin = ymint - .2*(ymaxt - ymint);
595:   ymax = ymaxt + .2*(ymaxt - ymint);
596: #if defined(PETSC_HAVE_X11)
597:   PetscDrawOpenX(PETSC_COMM_WORLD,0,"meshes",PETSC_DECIDE,PETSC_DECIDE,700,700,&draw);
598:   PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);
599:   PetscDrawZoom(draw,DrawPatch,&zctx);
600:   VecRestoreArray(v,&va);
601:   PetscDrawDestroy(draw);
602: #endif
603:   return(0);
604: }

606: /* crude mappings from rectangular arrays to the true geometry. These are ONLY for testing!
607:    they will not be used the actual code */
608: PetscErrorCode FAMapRegion3(FA fa,Vec g)
609: {
611:   PetscReal      R = 1.0,Rscale,Ascale;
612:   PetscInt       i,k,x,y,m,n;
613:   Field          **ga;

616:   Rscale = R/(fa->r2-1);
617:   Ascale = 2.0*PETSC_PI/(3.0*(fa->p1 - fa->p2 - 1));

619:   FAGetGlobalCorners(fa,2,&x,&y,&m,&n);
620:   FAGetGlobalArray(fa,g,2,&ga);
621:   for (k=y; k<y+n; k++) {
622:     for (i=x; i<x+m; i++) {
623:       ga[k][i].X = (R + k*Rscale)*PetscCosScalar(1.*PETSC_PI/6. + i*Ascale);
624:       ga[k][i].Y = (R + k*Rscale)*PetscSinScalar(1.*PETSC_PI/6. + i*Ascale) - 4.*R;
625:     }
626:   }
627:   FARestoreGlobalArray(fa,g,2,&ga);
628:   return(0);
629: }

631: PetscErrorCode FAMapRegion2(FA fa,Vec g)
632: {
634:   PetscReal      R = 1.0,Rscale,Ascale;
635:   PetscInt       i,k,x,y,m,n;
636:   Field          **ga;

639:   Rscale = R/(fa->r2-1);
640:   Ascale = 2.0*PETSC_PI/fa->p2;

642:   FAGetGlobalCorners(fa,1,&x,&y,&m,&n);
643:   FAGetGlobalArray(fa,g,1,&ga);
644:   for (k=y; k<y+n; k++) {
645:     for (i=x; i<x+m; i++) {
646:       ga[k][i].X = (R + k*Rscale)*PetscCosScalar(i*Ascale - PETSC_PI/2.0);
647:       ga[k][i].Y = (R + k*Rscale)*PetscSinScalar(i*Ascale - PETSC_PI/2.0);
648:     }
649:   }
650:   FARestoreGlobalArray(fa,g,1,&ga);
651:   return(0);
652: }

654: PetscErrorCode FAMapRegion1(FA fa,Vec g)
655: {
657:   PetscReal      R = 1.0,Rscale,Ascale1,Ascale3;
658:   PetscInt       i,k,x,y,m,n;
659:   Field          **ga;

662:   Rscale  = R/(fa->r1-1);
663:   Ascale1 = 2.0*PETSC_PI/fa->p2;
664:   Ascale3 = 2.0*PETSC_PI/(3.0*(fa->p1 - fa->p2 - 1));

666:   FAGetGlobalCorners(fa,0,&x,&y,&m,&n);
667:   FAGetGlobalArray(fa,g,0,&ga);

669:   /* This mapping is WRONG! Not sure how to do it so I've done a poor job of
670:      it. You can see that the grid connections are correct. */
671:   for (k=y; k<y+n; k++) {
672:     for (i=x; i<x+m; i++) {
673:       if (i < (fa->p1-fa->p2)/2) {
674:         ga[k][i].X = (2.0*R + k*Rscale)*PetscCosScalar(i*Ascale3);
675:         ga[k][i].Y = (2.0*R + k*Rscale)*PetscSinScalar(i*Ascale3) - 4.*R;
676:       } else if (i > fa->p2 + (fa->p1 - fa->p2)/2) {
677:         ga[k][i].X = (2.0*R + k*Rscale)*PetscCosScalar(PETSC_PI+i*Ascale3);
678:         ga[k][i].Y = (2.0*R + k*Rscale)*PetscSinScalar(PETSC_PI+i*Ascale3) - 4.*R;
679:       } else {
680:         ga[k][i].X = (2.*R + k*Rscale)*PetscCosScalar((i-(fa->p1-fa->p2)/2)*Ascale1 - PETSC_PI/2.0);
681:         ga[k][i].Y = (2.*R + k*Rscale)*PetscSinScalar((i-(fa->p1-fa->p2)/2)*Ascale1 - PETSC_PI/2.0);
682:       }
683:     }
684:   }
685:   FARestoreGlobalArray(fa,g,0,&ga);
686:   return(0);
687: }

689: /* Simple test to check that the ghost points are properly updated */
690: PetscErrorCode FATest(FA fa)
691: {
693:   Vec            l,g;
694:   Field          **la;
695:   PetscInt       x,y,m,n,j,i,k,p;
696:   PetscMPIInt    rank;

699:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

701:   FAGetGlobalVector(fa,&g);
702:   FAGetLocalVector(fa,&l);

704:   /* fill up global vector of one region at a time with ITS logical coordinates, then update LOCAL
705:      vector; print local vectors to confirm they are correctly filled */
706:   for (j=0; j<3; j++) {
707:     VecSet(g,0.0);
708:     FAGetGlobalCorners(fa,j,&x,&y,&m,&n);
709:     PetscPrintf(PETSC_COMM_WORLD,"\nFilling global region %d, showing local results \n",j+1);
710:     FAGetGlobalArray(fa,g,j,&la);
711:     for (k=y; k<y+n; k++) {
712:       for (i=x; i<x+m; i++) {
713:         la[k][i].X = i;
714:         la[k][i].Y = k;
715:       }
716:     }
717:     FARestoreGlobalArray(fa,g,j,&la);

719:     FAGlobalToLocal(fa,g,l);
720:     DrawFA(fa,g);
721:     DrawFA(fa,l);

723:     for (p=0; p<3; p++) {
724:       FAGetLocalCorners(fa,p,&x,&y,&m,&n);
725:       FAGetLocalArray(fa,l,p,&la);
726:       PetscSynchronizedPrintf(PETSC_COMM_WORLD,"\n[%d] Local array for region %d \n",rank,p+1);
727:       for (k=y+n-1; k>=y; k--) { /* print in reverse order to match diagram in paper */
728:         for (i=x; i<x+m; i++) {
729:           PetscSynchronizedPrintf(PETSC_COMM_WORLD,"(%G,%G) ",la[k][i].X,la[k][i].Y);
730:         }
731:         PetscSynchronizedPrintf(PETSC_COMM_WORLD,"\n");
732:       }
733:       FARestoreLocalArray(fa,l,p,&la);
734:       PetscSynchronizedFlush(PETSC_COMM_WORLD);
735:     }
736:   }
737:   VecDestroy(g);
738:   VecDestroy(l);
739:   return(0);
740: }

744: int main(int argc,char **argv)
745: {
746:   FA             fa;
748:   Vec            g,l;

750:   PetscInitialize(&argc,&argv,0,help);

752:   FACreate(&fa);
753:   /* FATest(fa);*/

755:   FAGetGlobalVector(fa,&g);
756:   FAGetLocalVector(fa,&l);

758:   FAMapRegion1(fa,g);
759:   FAMapRegion2(fa,g);
760:   FAMapRegion3(fa,g);

762:   FAGlobalToLocal(fa,g,l);
763:   DrawFA(fa,g);
764:   DrawFA(fa,l);

766:   VecDestroy(g);
767:   VecDestroy(l);
768:   FADestroy(fa);

770:   PetscFinalize();
771:   return 0;
772: }