Actual source code: da2.c

petsc-3.3-p7 2013-05-11
  2: #include <petsc-private/daimpl.h>    /*I   "petscdmda.h"   I*/

  6: PetscErrorCode DMView_DA_2d(DM da,PetscViewer viewer)
  7: {
  9:   PetscMPIInt    rank;
 10:   PetscBool      iascii,isdraw,isbinary;
 11:   DM_DA          *dd = (DM_DA*)da->data;
 12: #if defined(PETSC_HAVE_MATLAB_ENGINE)
 13:   PetscBool      ismatlab;
 14: #endif

 17:   MPI_Comm_rank(((PetscObject)da)->comm,&rank);

 19:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
 20:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
 21:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
 22: #if defined(PETSC_HAVE_MATLAB_ENGINE)
 23:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATLAB,&ismatlab);
 24: #endif
 25:   if (iascii) {
 26:     PetscViewerFormat format;

 28:     PetscViewerGetFormat(viewer, &format);
 29:     if (format != PETSC_VIEWER_ASCII_VTK && format != PETSC_VIEWER_ASCII_VTK_CELL) {
 30:       DMDALocalInfo info;
 31:       DMDAGetLocalInfo(da,&info);
 32:       PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
 33:       PetscViewerASCIISynchronizedPrintf(viewer,"Processor [%d] M %D N %D m %D n %D w %D s %D\n",rank,dd->M,dd->N,dd->m,dd->n,dd->w,dd->s);
 34:       PetscViewerASCIISynchronizedPrintf(viewer,"X range of indices: %D %D, Y range of indices: %D %D\n",info.xs,info.xs+info.xm,info.ys,info.ys+info.ym);
 35:       PetscViewerFlush(viewer);
 36:       PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
 37:     } else {
 38:       DMView_DA_VTK(da,viewer);
 39:     }
 40:   } else if (isdraw) {
 41:     PetscDraw  draw;
 42:     double     ymin = -1*dd->s-1,ymax = dd->N+dd->s;
 43:     double     xmin = -1*dd->s-1,xmax = dd->M+dd->s;
 44:     double     x,y;
 45:     PetscInt   base,*idx;
 46:     char       node[10];
 47:     PetscBool  isnull;
 48: 
 49:     PetscViewerDrawGetDraw(viewer,0,&draw);
 50:     PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
 51:     if (!dd->coordinates) {
 52:       PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);
 53:     }
 54:     PetscDrawSynchronizedClear(draw);

 56:     /* first processor draw all node lines */
 57:     if (!rank) {
 58:       ymin = 0.0; ymax = dd->N - 1;
 59:       for (xmin=0; xmin<dd->M; xmin++) {
 60:         PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_BLACK);
 61:       }
 62:       xmin = 0.0; xmax = dd->M - 1;
 63:       for (ymin=0; ymin<dd->N; ymin++) {
 64:         PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_BLACK);
 65:       }
 66:     }
 67:     PetscDrawSynchronizedFlush(draw);
 68:     PetscDrawPause(draw);

 70:     /* draw my box */
 71:     ymin = dd->ys; ymax = dd->ye - 1; xmin = dd->xs/dd->w;
 72:     xmax =(dd->xe-1)/dd->w;
 73:     PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_RED);
 74:     PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_RED);
 75:     PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_RED);
 76:     PetscDrawLine(draw,xmax,ymin,xmax,ymax,PETSC_DRAW_RED);

 78:     /* put in numbers */
 79:     base = (dd->base)/dd->w;
 80:     for (y=ymin; y<=ymax; y++) {
 81:       for (x=xmin; x<=xmax; x++) {
 82:         sprintf(node,"%d",(int)base++);
 83:         PetscDrawString(draw,x,y,PETSC_DRAW_BLACK,node);
 84:       }
 85:     }

 87:     PetscDrawSynchronizedFlush(draw);
 88:     PetscDrawPause(draw);
 89:     /* overlay ghost numbers, useful for error checking */
 90:     /* put in numbers */

 92:     base = 0; idx = dd->idx;
 93:     ymin = dd->Ys; ymax = dd->Ye; xmin = dd->Xs; xmax = dd->Xe;
 94:     for (y=ymin; y<ymax; y++) {
 95:       for (x=xmin; x<xmax; x++) {
 96:         if ((base % dd->w) == 0) {
 97:           sprintf(node,"%d",(int)(idx[base]/dd->w));
 98:           PetscDrawString(draw,x/dd->w,y,PETSC_DRAW_BLUE,node);
 99:         }
100:         base++;
101:       }
102:     }
103:     PetscDrawSynchronizedFlush(draw);
104:     PetscDrawPause(draw);
105:   } else if (isbinary){
106:     DMView_DA_Binary(da,viewer);
107: #if defined(PETSC_HAVE_MATLAB_ENGINE)
108:   } else if (ismatlab) {
109:     DMView_DA_Matlab(da,viewer);
110: #endif
111:   } else SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for DMDA 1d",((PetscObject)viewer)->type_name);
112:   return(0);
113: }

115: /*
116:       M is number of grid points 
117:       m is number of processors

119: */
122: PetscErrorCode  DMDASplitComm2d(MPI_Comm comm,PetscInt M,PetscInt N,PetscInt sw,MPI_Comm *outcomm)
123: {
125:   PetscInt       m,n = 0,x = 0,y = 0;
126:   PetscMPIInt    size,csize,rank;

129:   MPI_Comm_size(comm,&size);
130:   MPI_Comm_rank(comm,&rank);

132:   csize = 4*size;
133:   do {
134:     if (csize % 4) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cannot split communicator of size %d tried %d %D %D",size,csize,x,y);
135:     csize   = csize/4;
136: 
137:     m = (PetscInt)(0.5 + sqrt(((double)M)*((double)csize)/((double)N)));
138:     if (!m) m = 1;
139:     while (m > 0) {
140:       n = csize/m;
141:       if (m*n == csize) break;
142:       m--;
143:     }
144:     if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;}

146:     x = M/m + ((M % m) > ((csize-1) % m));
147:     y = (N + (csize-1)/m)/n;
148:   } while ((x < 4 || y < 4) && csize > 1);
149:   if (size != csize) {
150:     MPI_Group    entire_group,sub_group;
151:     PetscMPIInt  i,*groupies;

153:     MPI_Comm_group(comm,&entire_group);
154:     PetscMalloc(csize*sizeof(PetscInt),&groupies);
155:     for (i=0; i<csize; i++) {
156:       groupies[i] = (rank/csize)*csize + i;
157:     }
158:     MPI_Group_incl(entire_group,csize,groupies,&sub_group);
159:     PetscFree(groupies);
160:     MPI_Comm_create(comm,sub_group,outcomm);
161:     MPI_Group_free(&entire_group);
162:     MPI_Group_free(&sub_group);
163:     PetscInfo1(0,"DMDASplitComm2d:Creating redundant coarse problems of size %d\n",csize);
164:   } else {
165:     *outcomm = comm;
166:   }
167:   return(0);
168: }

172: static PetscErrorCode DMDAFunction(DM dm,Vec x,Vec F)
173: {
175:   Vec            localX;
176: 
178:   DMGetLocalVector(dm,&localX);
179:   DMGlobalToLocalBegin(dm,x,INSERT_VALUES,localX);
180:   DMGlobalToLocalEnd(dm,x,INSERT_VALUES,localX);
181:   DMDAComputeFunction1(dm,localX,F,dm->ctx);
182:   DMRestoreLocalVector(dm,&localX);
183:   return(0);
184: }

188: /*@C
189:        DMDASetLocalFunction - Caches in a DM a local function. 

191:    Logically Collective on DMDA

193:    Input Parameter:
194: +  da - initial distributed array
195: -  lf - the local function

197:    Level: intermediate

199:    Notes: 
200:       If you use SNESSetDM() and did not use SNESSetFunction() then the context passed to your function is the context set with DMSetApplicationContext()

202:    Developer Notes: It is possibly confusing which context is passed to the user function, it would be nice to unify them somehow.

204: .keywords:  distributed array, refine

206: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunctioni()
207: @*/
208: PetscErrorCode  DMDASetLocalFunction(DM da,DMDALocalFunction1 lf)
209: {
211:   DM_DA          *dd = (DM_DA*)da->data;

215:   DMSetFunction(da,DMDAFunction);
216:   dd->lf       = lf;
217:   return(0);
218: }

222: /*@C
223:        DMDASetLocalFunctioni - Caches in a DM a local function that evaluates a single component

225:    Logically Collective on DMDA

227:    Input Parameter:
228: +  da - initial distributed array
229: -  lfi - the local function

231:    Level: intermediate

233: .keywords:  distributed array, refine

235: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction()
236: @*/
237: PetscErrorCode  DMDASetLocalFunctioni(DM da,PetscErrorCode (*lfi)(DMDALocalInfo*,MatStencil*,void*,PetscScalar*,void*))
238: {
239:   DM_DA          *dd = (DM_DA*)da->data;
242:   dd->lfi = lfi;
243:   return(0);
244: }

248: /*@C
249:        DMDASetLocalFunctionib - Caches in a DM a block local function that evaluates a single component

251:    Logically Collective on DMDA

253:    Input Parameter:
254: +  da - initial distributed array
255: -  lfi - the local function

257:    Level: intermediate

259: .keywords:  distributed array, refine

261: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction()
262: @*/
263: PetscErrorCode  DMDASetLocalFunctionib(DM da,PetscErrorCode (*lfi)(DMDALocalInfo*,MatStencil*,void*,PetscScalar*,void*))
264: {
265:   DM_DA          *dd = (DM_DA*)da->data;
268:   dd->lfib = lfi;
269:   return(0);
270: }

274: PetscErrorCode DMDASetLocalAdicFunction_Private(DM da,DMDALocalFunction1 ad_lf)
275: {
276:   DM_DA          *dd = (DM_DA*)da->data;
279:   dd->adic_lf = ad_lf;
280:   return(0);
281: }

283: /*MC
284:        DMDASetLocalAdicFunctioni - Caches in a DM a local functioni computed by ADIC/ADIFOR

286:    Synopsis:
287:    PetscErrorCode DMDASetLocalAdicFunctioni(DM da,PetscInt (ad_lf*)(DMDALocalInfo*,MatStencil*,void*,void*,void*)
288:    
289:    Logically Collective on DMDA

291:    Input Parameter:
292: +  da - initial distributed array
293: -  ad_lfi - the local function as computed by ADIC/ADIFOR

295:    Level: intermediate

297: .keywords:  distributed array, refine

299: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction(),
300:           DMDASetLocalJacobian(), DMDASetLocalFunctioni()
301: M*/

305: PetscErrorCode DMDASetLocalAdicFunctioni_Private(DM da,PetscErrorCode (*ad_lfi)(DMDALocalInfo*,MatStencil*,void*,void*,void*))
306: {
307:   DM_DA          *dd = (DM_DA*)da->data;
310:   dd->adic_lfi = ad_lfi;
311:   return(0);
312: }

314: /*MC
315:        DMDASetLocalAdicMFFunctioni - Caches in a DM a local functioni computed by ADIC/ADIFOR

317:    Synopsis:
318:    PetscErrorCode  DMDASetLocalAdicFunctioni(DM da,int (ad_lf*)(DMDALocalInfo*,MatStencil*,void*,void*,void*)
319:    
320:    Logically Collective on DMDA

322:    Input Parameter:
323: +  da - initial distributed array
324: -  admf_lfi - the local matrix-free function as computed by ADIC/ADIFOR

326:    Level: intermediate

328: .keywords:  distributed array, refine

330: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction(),
331:           DMDASetLocalJacobian(), DMDASetLocalFunctioni()
332: M*/

336: PetscErrorCode DMDASetLocalAdicMFFunctioni_Private(DM da,PetscErrorCode (*admf_lfi)(DMDALocalInfo*,MatStencil*,void*,void*,void*))
337: {
338:   DM_DA          *dd = (DM_DA*)da->data;
341:   dd->adicmf_lfi = admf_lfi;
342:   return(0);
343: }

345: /*MC
346:        DMDASetLocalAdicFunctionib - Caches in a DM a block local functioni computed by ADIC/ADIFOR

348:    Synopsis:
349:    PetscErrorCode DMDASetLocalAdicFunctionib(DM da,PetscInt (ad_lf*)(DMDALocalInfo*,MatStencil*,void*,void*,void*)
350:    
351:    Logically Collective on DMDA

353:    Input Parameter:
354: +  da - initial distributed array
355: -  ad_lfi - the local function as computed by ADIC/ADIFOR

357:    Level: intermediate

359: .keywords:  distributed array, refine

361: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction(),
362:           DMDASetLocalJacobian(), DMDASetLocalFunctionib()
363: M*/

367: PetscErrorCode DMDASetLocalAdicFunctionib_Private(DM da,PetscErrorCode (*ad_lfi)(DMDALocalInfo*,MatStencil*,void*,void*,void*))
368: {
369:   DM_DA          *dd = (DM_DA*)da->data;
372:   dd->adic_lfib = ad_lfi;
373:   return(0);
374: }

376: /*MC
377:        DMDASetLocalAdicMFFunctionib - Caches in a DM a block local functioni computed by ADIC/ADIFOR

379:    Synopsis:
380:    PetscErrorCode  DMDASetLocalAdicFunctionib(DM da,int (ad_lf*)(DMDALocalInfo*,MatStencil*,void*,void*,void*)

382:    Logically Collective on DMDA

384:    Input Parameter:
385: +  da - initial distributed array
386: -  admf_lfi - the local matrix-free function as computed by ADIC/ADIFOR

388:    Level: intermediate

390: .keywords:  distributed array, refine

392: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction(),
393:           DMDASetLocalJacobian(), DMDASetLocalFunctionib()
394: M*/

398: PetscErrorCode DMDASetLocalAdicMFFunctionib_Private(DM da,PetscErrorCode (*admf_lfi)(DMDALocalInfo*,MatStencil*,void*,void*,void*))
399: {
400:   DM_DA          *dd = (DM_DA*)da->data;
403:   dd->adicmf_lfib = admf_lfi;
404:   return(0);
405: }

407: /*MC
408:        DMDASetLocalAdicMFFunction - Caches in a DM a local function computed by ADIC/ADIFOR

410:    Synopsis:
411:    PetscErrorCode DMDASetLocalAdicMFFunction(DM da,DMDALocalFunction1 ad_lf)

413:    Logically Collective on DMDA

415:    Input Parameter:
416: +  da - initial distributed array
417: -  ad_lf - the local function as computed by ADIC/ADIFOR

419:    Level: intermediate

421: .keywords:  distributed array, refine

423: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction(),
424:           DMDASetLocalJacobian()
425: M*/

429: PetscErrorCode DMDASetLocalAdicMFFunction_Private(DM da,DMDALocalFunction1 ad_lf)
430: {
431:   DM_DA          *dd = (DM_DA*)da->data;
434:   dd->adicmf_lf = ad_lf;
435:   return(0);
436: }

440: PetscErrorCode DMDAJacobianLocal(DM dm,Vec x,Mat A,Mat B, MatStructure *str)
441: {
443:   Vec            localX;
444: 
446:   DMGetLocalVector(dm,&localX);
447:   DMGlobalToLocalBegin(dm,x,INSERT_VALUES,localX);
448:   DMGlobalToLocalEnd(dm,x,INSERT_VALUES,localX);
449:   MatFDColoringApply(B,dm->fd,localX,str,dm);
450:   DMRestoreLocalVector(dm,&localX);
451:   /* Assemble true Jacobian; if it is different */
452:   if (A != B) {
453:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
454:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
455:   }
456:   MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
457:   *str = SAME_NONZERO_PATTERN;
458:   return(0);
459: }


464: static PetscErrorCode DMDAJacobian(DM dm,Vec x,Mat A,Mat B, MatStructure *str)
465: {
467:   Vec            localX;
468: 
470:   DMGetLocalVector(dm,&localX);
471:   DMGlobalToLocalBegin(dm,x,INSERT_VALUES,localX);
472:   DMGlobalToLocalEnd(dm,x,INSERT_VALUES,localX);
473:   DMDAComputeJacobian1(dm,localX,B,dm->ctx);
474:   DMRestoreLocalVector(dm,&localX);
475:   /* Assemble true Jacobian; if it is different */
476:   if (A != B) {
477:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
478:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
479:   }
480:   MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
481:   *str = SAME_NONZERO_PATTERN;
482:   return(0);
483: }

485: /*@C
486:        DMDASetLocalJacobian - Caches in a DM a local Jacobian computation function

488:    Logically Collective on DMDA

490:    
491:    Input Parameter:
492: +  da - initial distributed array
493: -  lj - the local Jacobian

495:    Level: intermediate

497: .keywords:  distributed array, refine

499: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction()
500: @*/
503: PetscErrorCode  DMDASetLocalJacobian(DM da,DMDALocalFunction1 lj)
504: {
506:   DM_DA          *dd = (DM_DA*)da->data;

510:   DMSetJacobian(da,DMDAJacobian);
511:   dd->lj    = lj;
512:   return(0);
513: }

517: /*@C
518:        DMDAGetLocalFunction - Gets from a DM a local function and its ADIC/ADIFOR Jacobian

520:    Note Collective

522:    Input Parameter:
523: .  da - initial distributed array

525:    Output Parameter:
526: .  lf - the local function

528:    Level: intermediate

530: .keywords:  distributed array, refine

532: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalJacobian(), DMDASetLocalFunction()
533: @*/
534: PetscErrorCode  DMDAGetLocalFunction(DM da,DMDALocalFunction1 *lf)
535: {
536:   DM_DA *dd = (DM_DA*)da->data;
539:   if (lf) *lf = dd->lf;
540:   return(0);
541: }

545: /*@C
546:        DMDAGetLocalJacobian - Gets from a DM a local jacobian

548:    Not Collective

550:    Input Parameter:
551: .  da - initial distributed array

553:    Output Parameter:
554: .  lj - the local jacobian

556:    Level: intermediate

558: .keywords:  distributed array, refine

560: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalJacobian()
561: @*/
562: PetscErrorCode  DMDAGetLocalJacobian(DM da,DMDALocalFunction1 *lj)
563: {
564:   DM_DA *dd = (DM_DA*)da->data;
567:   if (lj) *lj = dd->lj;
568:   return(0);
569: }

573: /*@
574:     DMDAComputeFunction - Evaluates a user provided function on each processor that 
575:         share a DMDA

577:    Input Parameters:
578: +    da - the DM that defines the grid
579: .    vu - input vector
580: .    vfu - output vector 
581: -    w - any user data

583:     Notes: Does NOT do ghost updates on vu upon entry

585:            This should eventually replace DMDAComputeFunction1

587:     Level: advanced

589: .seealso: DMDAComputeJacobian1WithAdic()

591: @*/
592: PetscErrorCode  DMDAComputeFunction(DM da,PetscErrorCode (*lf)(void),Vec vu,Vec vfu,void *w)
593: {
595:   void           *u,*fu;
596:   DMDALocalInfo  info;
597:   PetscErrorCode (*f)(DMDALocalInfo*,void*,void*,void*) = (PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))lf;
598: 
600:   DMDAGetLocalInfo(da,&info);
601:   DMDAVecGetArray(da,vu,&u);
602:   DMDAVecGetArray(da,vfu,&fu);

604:   (*f)(&info,u,fu,w);

606:   DMDAVecRestoreArray(da,vu,&u);
607:   DMDAVecRestoreArray(da,vfu,&fu);
608:   return(0);
609: }

613: /*@C 
614:    DMDAComputeFunctionLocal - This is a universal function evaluation routine for
615:    a local DM function.

617:    Collective on DMDA

619:    Input Parameters:
620: +  da - the DM context
621: .  func - The local function
622: .  X - input vector
623: .  F - function vector
624: -  ctx - A user context

626:    Level: intermediate

628: .seealso: DMDASetLocalFunction(), DMDASetLocalJacobian(), DMDASetLocalAdicFunction(), DMDASetLocalAdicMFFunction(),
629:           SNESSetFunction(), SNESSetJacobian()

631: @*/
632: PetscErrorCode  DMDAComputeFunctionLocal(DM da, DMDALocalFunction1 func, Vec X, Vec F, void *ctx)
633: {
634:   Vec            localX;
635:   DMDALocalInfo  info;
636:   void           *u;
637:   void           *fu;

641:   DMGetLocalVector(da,&localX);
642:   /*
643:      Scatter ghost points to local vector, using the 2-step process
644:         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
645:   */
646:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
647:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
648:   DMDAGetLocalInfo(da,&info);
649:   DMDAVecGetArray(da,localX,&u);
650:   DMDAVecGetArray(da,F,&fu);
651:   (*func)(&info,u,fu,ctx);
652:   DMDAVecRestoreArray(da,localX,&u);
653:   DMDAVecRestoreArray(da,F,&fu);
654:   DMRestoreLocalVector(da,&localX);
655:   return(0);
656: }

660: /*@C 
661:    DMDAComputeFunctionLocalGhost - This is a universal function evaluation routine for
662:    a local DM function, but the ghost values of the output are communicated and added.

664:    Collective on DMDA

666:    Input Parameters:
667: +  da - the DM context
668: .  func - The local function
669: .  X - input vector
670: .  F - function vector
671: -  ctx - A user context

673:    Level: intermediate

675: .seealso: DMDASetLocalFunction(), DMDASetLocalJacobian(), DMDASetLocalAdicFunction(), DMDASetLocalAdicMFFunction(),
676:           SNESSetFunction(), SNESSetJacobian()

678: @*/
679: PetscErrorCode  DMDAComputeFunctionLocalGhost(DM da, DMDALocalFunction1 func, Vec X, Vec F, void *ctx)
680: {
681:   Vec            localX, localF;
682:   DMDALocalInfo  info;
683:   void           *u;
684:   void           *fu;

688:   DMGetLocalVector(da,&localX);
689:   DMGetLocalVector(da,&localF);
690:   /*
691:      Scatter ghost points to local vector, using the 2-step process
692:         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
693:   */
694:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
695:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
696:   VecSet(F, 0.0);
697:   VecSet(localF, 0.0);
698:   DMDAGetLocalInfo(da,&info);
699:   DMDAVecGetArray(da,localX,&u);
700:   DMDAVecGetArray(da,localF,&fu);
701:   (*func)(&info,u,fu,ctx);
702:   DMLocalToGlobalBegin(da,localF,ADD_VALUES,F);
703:   DMLocalToGlobalEnd(da,localF,ADD_VALUES,F);
704:   DMDAVecRestoreArray(da,localX,&u);
705:   DMDAVecRestoreArray(da,localF,&fu);
706:   DMRestoreLocalVector(da,&localX);
707:   DMRestoreLocalVector(da,&localF);
708:   return(0);
709: }

713: /*@
714:     DMDAComputeFunction1 - Evaluates a user provided function on each processor that 
715:         share a DMDA

717:    Input Parameters:
718: +    da - the DM that defines the grid
719: .    vu - input vector
720: .    vfu - output vector 
721: -    w - any user data

723:     Notes: Does NOT do ghost updates on vu upon entry

725:     Level: advanced

727: .seealso: DMDAComputeJacobian1WithAdic()

729: @*/
730: PetscErrorCode  DMDAComputeFunction1(DM da,Vec vu,Vec vfu,void *w)
731: {
733:   void           *u,*fu;
734:   DMDALocalInfo  info;
735:   DM_DA          *dd = (DM_DA*)da->data;
736: 
738:   if (!dd->lf) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_NULL,"DMDASetLocalFunction() never called");
739:   DMDAGetLocalInfo(da,&info);
740:   DMDAVecGetArray(da,vu,&u);
741:   DMDAVecGetArray(da,vfu,&fu);

743:   CHKMEMQ;
744:   (*dd->lf)(&info,u,fu,w);
745:   CHKMEMQ;

747:   DMDAVecRestoreArray(da,vu,&u);
748:   DMDAVecRestoreArray(da,vfu,&fu);
749:   return(0);
750: }

754: PetscErrorCode  DMDAComputeFunctioniTest1(DM da,void *w)
755: {
756:   Vec            vu,fu,fui;
758:   PetscInt       i,n;
759:   PetscScalar    *ui;
760:   PetscRandom    rnd;
761:   PetscReal      norm;

764:   DMGetLocalVector(da,&vu);
765:   PetscRandomCreate(PETSC_COMM_SELF,&rnd);
766:   PetscRandomSetFromOptions(rnd);
767:   VecSetRandom(vu,rnd);
768:   PetscRandomDestroy(&rnd);

770:   DMGetGlobalVector(da,&fu);
771:   DMGetGlobalVector(da,&fui);
772: 
773:   DMDAComputeFunction1(da,vu,fu,w);

775:   VecGetArray(fui,&ui);
776:   VecGetLocalSize(fui,&n);
777:   for (i=0; i<n; i++) {
778:     DMDAComputeFunctioni1(da,i,vu,ui+i,w);
779:   }
780:   VecRestoreArray(fui,&ui);

782:   VecAXPY(fui,-1.0,fu);
783:   VecNorm(fui,NORM_2,&norm);
784:   PetscPrintf(((PetscObject)da)->comm,"Norm of difference in vectors %G\n",norm);
785:   VecView(fu,0);
786:   VecView(fui,0);

788:   DMRestoreLocalVector(da,&vu);
789:   DMRestoreGlobalVector(da,&fu);
790:   DMRestoreGlobalVector(da,&fui);
791:   return(0);
792: }

796: /*@
797:     DMDAComputeFunctioni1 - Evaluates a user provided point-wise function

799:    Input Parameters:
800: +    da - the DM that defines the grid
801: .    i - the component of the function we wish to compute (must be local)
802: .    vu - input vector
803: .    vfu - output value
804: -    w - any user data

806:     Notes: Does NOT do ghost updates on vu upon entry

808:     Level: advanced

810: .seealso: DMDAComputeJacobian1WithAdic()

812: @*/
813: PetscErrorCode  DMDAComputeFunctioni1(DM da,PetscInt i,Vec vu,PetscScalar *vfu,void *w)
814: {
816:   void           *u;
817:   DMDALocalInfo  info;
818:   MatStencil     stencil;
819:   DM_DA          *dd = (DM_DA*)da->data;
820: 

823:   DMDAGetLocalInfo(da,&info);
824:   DMDAVecGetArray(da,vu,&u);

826:   /* figure out stencil value from i */
827:   stencil.c = i % info.dof;
828:   stencil.i = (i % (info.xm*info.dof))/info.dof;
829:   stencil.j = (i % (info.xm*info.ym*info.dof))/(info.xm*info.dof);
830:   stencil.k = i/(info.xm*info.ym*info.dof);

832:   (*dd->lfi)(&info,&stencil,u,vfu,w);

834:   DMDAVecRestoreArray(da,vu,&u);
835:   return(0);
836: }

840: /*@
841:     DMDAComputeFunctionib1 - Evaluates a user provided point-block function

843:    Input Parameters:
844: +    da - the DM that defines the grid
845: .    i - the component of the function we wish to compute (must be local)
846: .    vu - input vector
847: .    vfu - output value
848: -    w - any user data

850:     Notes: Does NOT do ghost updates on vu upon entry

852:     Level: advanced

854: .seealso: DMDAComputeJacobian1WithAdic()

856: @*/
857: PetscErrorCode  DMDAComputeFunctionib1(DM da,PetscInt i,Vec vu,PetscScalar *vfu,void *w)
858: {
860:   void           *u;
861:   DMDALocalInfo  info;
862:   MatStencil     stencil;
863:   DM_DA          *dd = (DM_DA*)da->data;
864: 
866:   DMDAGetLocalInfo(da,&info);
867:   DMDAVecGetArray(da,vu,&u);

869:   /* figure out stencil value from i */
870:   stencil.c = i % info.dof;
871:   if (stencil.c) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Point-block functions can only be called for the entire block");
872:   stencil.i = (i % (info.xm*info.dof))/info.dof;
873:   stencil.j = (i % (info.xm*info.ym*info.dof))/(info.xm*info.dof);
874:   stencil.k = i/(info.xm*info.ym*info.dof);

876:   (*dd->lfib)(&info,&stencil,u,vfu,w);

878:   DMDAVecRestoreArray(da,vu,&u);
879:   return(0);
880: }

882: #if defined(new)
885: /*
886:   DMDAGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix where local
887:     function lives on a DMDA

889:         y ~= (F(u + ha) - F(u))/h, 
890:   where F = nonlinear function, as set by SNESSetFunction()
891:         u = current iterate
892:         h = difference interval
893: */
894: PetscErrorCode DMDAGetDiagonal_MFFD(DM da,Vec U,Vec a)
895: {
896:   PetscScalar    h,*aa,*ww,v;
897:   PetscReal      epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
899:   PetscInt       gI,nI;
900:   MatStencil     stencil;
901:   DMDALocalInfo  info;
902: 
904:   (*ctx->func)(0,U,a,ctx->funcctx);
905:   (*ctx->funcisetbase)(U,ctx->funcctx);

907:   VecGetArray(U,&ww);
908:   VecGetArray(a,&aa);
909: 
910:   nI = 0;
911:     h  = ww[gI];
912:     if (h == 0.0) h = 1.0;
913: #if !defined(PETSC_USE_COMPLEX)
914:     if (h < umin && h >= 0.0)      h = umin;
915:     else if (h < 0.0 && h > -umin) h = -umin;
916: #else
917:     if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0)     h = umin;
918:     else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
919: #endif
920:     h     *= epsilon;
921: 
922:     ww[gI] += h;
923:     (*ctx->funci)(i,w,&v,ctx->funcctx);
924:     aa[nI]  = (v - aa[nI])/h;
925:     ww[gI] -= h;
926:     nI++;
927:   }
928:   VecRestoreArray(U,&ww);
929:   VecRestoreArray(a,&aa);
930:   return(0);
931: }
932: #endif

934: #if defined(PETSC_HAVE_ADIC)
935: EXTERN_C_BEGIN
936: #include <adic/ad_utils.h>
937: EXTERN_C_END

941: /*@C
942:     DMDAComputeJacobian1WithAdic - Evaluates a adiC provided Jacobian function on each processor that 
943:         share a DMDA

945:    Input Parameters:
946: +    da - the DM that defines the grid
947: .    vu - input vector (ghosted)
948: .    J - output matrix
949: -    w - any user data

951:    Level: advanced

953:     Notes: Does NOT do ghost updates on vu upon entry

955: .seealso: DMDAComputeFunction1()

957: @*/
958: PetscErrorCode  DMDAComputeJacobian1WithAdic(DM da,Vec vu,Mat J,void *w)
959: {
961:   PetscInt       gtdof,tdof;
962:   PetscScalar    *ustart;
963:   DMDALocalInfo  info;
964:   void           *ad_u,*ad_f,*ad_ustart,*ad_fstart;
965:   ISColoring     iscoloring;

968:   DMDAGetLocalInfo(da,&info);

970:   PetscADResetIndep();

972:   /* get space for derivative objects.  */
973:   DMDAGetAdicArray(da,PETSC_TRUE,&ad_u,&ad_ustart,&gtdof);
974:   DMDAGetAdicArray(da,PETSC_FALSE,&ad_f,&ad_fstart,&tdof);
975:   VecGetArray(vu,&ustart);
976:   DMCreateColoring(da,IS_COLORING_GHOSTED,MATAIJ,&iscoloring);

978:   PetscADSetValueAndColor(ad_ustart,gtdof,iscoloring->colors,ustart);

980:   VecRestoreArray(vu,&ustart);
981:   ISColoringDestroy(&iscoloring);
982:   PetscADIncrementTotalGradSize(iscoloring->n);
983:   PetscADSetIndepDone();

985:   PetscLogEventBegin(DMDA_LocalADFunction,0,0,0,0);
986:   (*dd->adic_lf)(&info,ad_u,ad_f,w);
987:   PetscLogEventEnd(DMDA_LocalADFunction,0,0,0,0);

989:   /* stick the values into the matrix */
990:   MatSetValuesAdic(J,(PetscScalar**)ad_fstart);

992:   /* return space for derivative objects.  */
993:   DMDARestoreAdicArray(da,PETSC_TRUE,&ad_u,&ad_ustart,&gtdof);
994:   DMDARestoreAdicArray(da,PETSC_FALSE,&ad_f,&ad_fstart,&tdof);
995:   return(0);
996: }

1000: /*@C
1001:     DMDAMultiplyByJacobian1WithAdic - Applies an ADIC-provided Jacobian function to a vector on 
1002:     each processor that shares a DMDA.

1004:     Input Parameters:
1005: +   da - the DM that defines the grid
1006: .   vu - Jacobian is computed at this point (ghosted)
1007: .   v - product is done on this vector (ghosted)
1008: .   fu - output vector = J(vu)*v (not ghosted)
1009: -   w - any user data

1011:     Notes: 
1012:     This routine does NOT do ghost updates on vu upon entry.

1014:    Level: advanced

1016: .seealso: DMDAComputeFunction1()

1018: @*/
1019: PetscErrorCode  DMDAMultiplyByJacobian1WithAdic(DM da,Vec vu,Vec v,Vec f,void *w)
1020: {
1022:   PetscInt       i,gtdof,tdof;
1023:   PetscScalar    *avu,*av,*af,*ad_vustart,*ad_fstart;
1024:   DMDALocalInfo  info;
1025:   void           *ad_vu,*ad_f;

1028:   DMDAGetLocalInfo(da,&info);

1030:   /* get space for derivative objects.  */
1031:   DMDAGetAdicMFArray(da,PETSC_TRUE,&ad_vu,&ad_vustart,&gtdof);
1032:   DMDAGetAdicMFArray(da,PETSC_FALSE,&ad_f,&ad_fstart,&tdof);

1034:   /* copy input vector into derivative object */
1035:   VecGetArray(vu,&avu);
1036:   VecGetArray(v,&av);
1037:   for (i=0; i<gtdof; i++) {
1038:     ad_vustart[2*i]   = avu[i];
1039:     ad_vustart[2*i+1] = av[i];
1040:   }
1041:   VecRestoreArray(vu,&avu);
1042:   VecRestoreArray(v,&av);

1044:   PetscADResetIndep();
1045:   PetscADIncrementTotalGradSize(1);
1046:   PetscADSetIndepDone();

1048:   (*dd->adicmf_lf)(&info,ad_vu,ad_f,w);

1050:   /* stick the values into the vector */
1051:   VecGetArray(f,&af);
1052:   for (i=0; i<tdof; i++) {
1053:     af[i] = ad_fstart[2*i+1];
1054:   }
1055:   VecRestoreArray(f,&af);

1057:   /* return space for derivative objects.  */
1058:   DMDARestoreAdicMFArray(da,PETSC_TRUE,&ad_vu,&ad_vustart,&gtdof);
1059:   DMDARestoreAdicMFArray(da,PETSC_FALSE,&ad_f,&ad_fstart,&tdof);
1060:   return(0);
1061: }
1062: #endif

1066: /*@
1067:     DMDAComputeJacobian1 - Evaluates a local Jacobian function on each processor that 
1068:         share a DMDA

1070:    Input Parameters:
1071: +    da - the DM that defines the grid
1072: .    vu - input vector (ghosted)
1073: .    J - output matrix
1074: -    w - any user data

1076:     Notes: Does NOT do ghost updates on vu upon entry

1078:     Level: advanced

1080: .seealso: DMDAComputeFunction1()

1082: @*/
1083: PetscErrorCode  DMDAComputeJacobian1(DM da,Vec vu,Mat J,void *w)
1084: {
1086:   void           *u;
1087:   DMDALocalInfo  info;
1088:   DM_DA          *dd = (DM_DA*)da->data;

1091:   DMDAGetLocalInfo(da,&info);
1092:   DMDAVecGetArray(da,vu,&u);
1093:   (*dd->lj)(&info,u,J,w);
1094:   DMDAVecRestoreArray(da,vu,&u);
1095:   return(0);
1096: }


1101: /*
1102:     DMDAComputeJacobian1WithAdifor - Evaluates a ADIFOR provided Jacobian local function on each processor that 
1103:         share a DMDA

1105:    Input Parameters:
1106: +    da - the DM that defines the grid
1107: .    vu - input vector (ghosted)
1108: .    J - output matrix
1109: -    w - any user data

1111:     Notes: Does NOT do ghost updates on vu upon entry

1113: .seealso: DMDAComputeFunction1()

1115: */
1116: PetscErrorCode  DMDAComputeJacobian1WithAdifor(DM da,Vec vu,Mat J,void *w)
1117: {
1118:   PetscErrorCode  ierr;
1119:   PetscInt        i,Nc,N;
1120:   ISColoringValue *color;
1121:   DMDALocalInfo   info;
1122:   PetscScalar     *u,*g_u,*g_f,*f = 0,*p_u;
1123:   ISColoring      iscoloring;
1124:   DM_DA          *dd = (DM_DA*)da->data;
1125:   void            (*lf)(PetscInt*,DMDALocalInfo*,PetscScalar*,PetscScalar*,PetscInt*,PetscScalar*,PetscScalar*,PetscInt*,void*,PetscErrorCode*) =
1126:                   (void (*)(PetscInt*,DMDALocalInfo*,PetscScalar*,PetscScalar*,PetscInt*,PetscScalar*,PetscScalar*,PetscInt*,void*,PetscErrorCode*))*dd->adifor_lf;

1129:   DMCreateColoring(da,IS_COLORING_GHOSTED,MATAIJ,&iscoloring);
1130:   Nc   = iscoloring->n;
1131:   DMDAGetLocalInfo(da,&info);
1132:   N    = info.gxm*info.gym*info.gzm*info.dof;

1134:   /* get space for derivative objects.  */
1135:   PetscMalloc(Nc*info.gxm*info.gym*info.gzm*info.dof*sizeof(PetscScalar),&g_u);
1136:   PetscMemzero(g_u,Nc*info.gxm*info.gym*info.gzm*info.dof*sizeof(PetscScalar));
1137:   p_u   = g_u;
1138:   color = iscoloring->colors;
1139:   for (i=0; i<N; i++) {
1140:     p_u[*color++] = 1.0;
1141:     p_u          += Nc;
1142:   }
1143:   ISColoringDestroy(&iscoloring);
1144:   PetscMalloc2(Nc*info.xm*info.ym*info.zm*info.dof,PetscScalar,&g_f,info.xm*info.ym*info.zm*info.dof,PetscScalar,&f);

1146:   /* Seed the input array g_u with coloring information */
1147: 
1148:   VecGetArray(vu,&u);
1149:   (lf)(&Nc,&info,u,g_u,&Nc,f,g_f,&Nc,w,&ierr);
1150:   VecRestoreArray(vu,&u);

1152:   /* stick the values into the matrix */
1153:   /* PetscScalarView(Nc*info.xm*info.ym,g_f,0); */
1154:   MatSetValuesAdifor(J,Nc,g_f);

1156:   /* return space for derivative objects.  */
1157:   PetscFree(g_u);
1158:   PetscFree2(g_f,f);
1159:   return(0);
1160: }

1164: /*@C 
1165:    DMDAFormjacobianLocal - This is a universal Jacobian evaluation routine for
1166:    a local DM function.

1168:    Collective on DMDA

1170:    Input Parameters:
1171: +  da - the DM context
1172: .  func - The local function
1173: .  X - input vector
1174: .  J - Jacobian matrix
1175: -  ctx - A user context

1177:    Level: intermediate

1179: .seealso: DMDASetLocalFunction(), DMDASetLocalJacobian(), DMDASetLocalAdicFunction(), DMDASetLocalAdicMFFunction(),
1180:           SNESSetFunction(), SNESSetJacobian()

1182: @*/
1183: PetscErrorCode  DMDAFormJacobianLocal(DM da, DMDALocalFunction1 func, Vec X, Mat J, void *ctx)
1184: {
1185:   Vec            localX;
1186:   DMDALocalInfo  info;
1187:   void           *u;

1191:   DMGetLocalVector(da,&localX);
1192:   /*
1193:      Scatter ghost points to local vector, using the 2-step process
1194:         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
1195:   */
1196:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
1197:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
1198:   DMDAGetLocalInfo(da,&info);
1199:   DMDAVecGetArray(da,localX,&u);
1200:   (*func)(&info,u,J,ctx);
1201:   DMDAVecRestoreArray(da,localX,&u);
1202:   DMRestoreLocalVector(da,&localX);
1203:   return(0);
1204: }

1208: /*@C
1209:     DMDAMultiplyByJacobian1WithAD - Applies a Jacobian function supplied by ADIFOR or ADIC
1210:     to a vector on each processor that shares a DMDA.

1212:    Input Parameters:
1213: +    da - the DM that defines the grid
1214: .    vu - Jacobian is computed at this point (ghosted)
1215: .    v - product is done on this vector (ghosted)
1216: .    fu - output vector = J(vu)*v (not ghosted)
1217: -    w - any user data

1219:     Notes: 
1220:     This routine does NOT do ghost updates on vu and v upon entry.
1221:            
1222:     Automatically calls DMDAMultiplyByJacobian1WithAdifor() or DMDAMultiplyByJacobian1WithAdic()
1223:     depending on whether DMDASetLocalAdicMFFunction() or DMDASetLocalAdiforMFFunction() was called.

1225:    Level: advanced

1227: .seealso: DMDAComputeFunction1(), DMDAMultiplyByJacobian1WithAdifor(), DMDAMultiplyByJacobian1WithAdic()

1229: @*/
1230: PetscErrorCode  DMDAMultiplyByJacobian1WithAD(DM da,Vec u,Vec v,Vec f,void *w)
1231: {
1233:   DM_DA          *dd = (DM_DA*)da->data;

1236:   if (dd->adicmf_lf) {
1237: #if defined(PETSC_HAVE_ADIC)
1238:     DMDAMultiplyByJacobian1WithAdic(da,u,v,f,w);
1239: #else
1240:     SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP_SYS,"Requires ADIC to be installed and cannot use complex numbers");
1241: #endif
1242:   } else if (dd->adiformf_lf) {
1243:     DMDAMultiplyByJacobian1WithAdifor(da,u,v,f,w);
1244:   } else SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ORDER,"Must call DMDASetLocalAdiforMFFunction() or DMDASetLocalAdicMFFunction() before using");
1245:   return(0);
1246: }


1251: /*@C
1252:     DMDAMultiplyByJacobian1WithAdifor - Applies a ADIFOR provided Jacobian function on each processor that 
1253:         share a DM to a vector

1255:    Input Parameters:
1256: +    da - the DM that defines the grid
1257: .    vu - Jacobian is computed at this point (ghosted)
1258: .    v - product is done on this vector (ghosted)
1259: .    fu - output vector = J(vu)*v (not ghosted)
1260: -    w - any user data

1262:     Notes: Does NOT do ghost updates on vu and v upon entry

1264:    Level: advanced

1266: .seealso: DMDAComputeFunction1()

1268: @*/
1269: PetscErrorCode  DMDAMultiplyByJacobian1WithAdifor(DM da,Vec u,Vec v,Vec f,void *w)
1270: {
1272:   PetscScalar    *au,*av,*af,*awork;
1273:   Vec            work;
1274:   DMDALocalInfo  info;
1275:   DM_DA          *dd = (DM_DA*)da->data;
1276:   void           (*lf)(DMDALocalInfo*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,void*,PetscErrorCode*) =
1277:                  (void (*)(DMDALocalInfo*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,void*,PetscErrorCode*))*dd->adiformf_lf;

1280:   DMDAGetLocalInfo(da,&info);

1282:   DMGetGlobalVector(da,&work);
1283:   VecGetArray(u,&au);
1284:   VecGetArray(v,&av);
1285:   VecGetArray(f,&af);
1286:   VecGetArray(work,&awork);
1287:   (lf)(&info,au,av,awork,af,w,&ierr);
1288:   VecRestoreArray(u,&au);
1289:   VecRestoreArray(v,&av);
1290:   VecRestoreArray(f,&af);
1291:   VecRestoreArray(work,&awork);
1292:   DMRestoreGlobalVector(da,&work);

1294:   return(0);
1295: }

1299: PetscErrorCode  DMSetUp_DA_2D(DM da)
1300: {
1301:   DM_DA                  *dd = (DM_DA*)da->data;
1302:   const PetscInt         M            = dd->M;
1303:   const PetscInt         N            = dd->N;
1304:   PetscInt               m            = dd->m;
1305:   PetscInt               n            = dd->n;
1306:   const PetscInt         dof          = dd->w;
1307:   const PetscInt         s            = dd->s;
1308:   const DMDABoundaryType bx         = dd->bx;
1309:   const DMDABoundaryType by         = dd->by;
1310:   const DMDAStencilType  stencil_type = dd->stencil_type;
1311:   PetscInt               *lx           = dd->lx;
1312:   PetscInt               *ly           = dd->ly;
1313:   MPI_Comm               comm;
1314:   PetscMPIInt            rank,size;
1315:   PetscInt               xs,xe,ys,ye,x,y,Xs,Xe,Ys,Ye,start,end,IXs,IXe,IYs,IYe;
1316:   PetscInt               up,down,left,right,i,n0,n1,n2,n3,n5,n6,n7,n8,*idx,nn,*idx_cpy;
1317:   const PetscInt         *idx_full;
1318:   PetscInt               xbase,*bases,*ldims,j,x_t,y_t,s_t,base,count;
1319:   PetscInt               s_x,s_y; /* s proportionalized to w */
1320:   PetscInt               sn0 = 0,sn2 = 0,sn6 = 0,sn8 = 0;
1321:   Vec                    local,global;
1322:   VecScatter             ltog,gtol;
1323:   IS                     to,from,ltogis;
1324:   PetscErrorCode         ierr;

1327:   if (dof < 1) SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Must have 1 or more degrees of freedom per node: %D",dof);
1328:   if (s < 0) SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Stencil width cannot be negative: %D",s);
1329:   PetscObjectGetComm((PetscObject)da,&comm);
1330: #if !defined(PETSC_USE_64BIT_INDICES)
1331:   if (((Petsc64bitInt) M)*((Petsc64bitInt) N)*((Petsc64bitInt) dof) > (Petsc64bitInt) PETSC_MPI_INT_MAX) SETERRQ3(comm,PETSC_ERR_INT_OVERFLOW,"Mesh of %D by %D by %D (dof) is too large for 32 bit indices",M,N,dof);
1332: #endif

1334:   if (dof < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Must have 1 or more degrees of freedom per node: %D",dof);
1335:   if (s < 0) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Stencil width cannot be negative: %D",s);

1337:   MPI_Comm_size(comm,&size);
1338:   MPI_Comm_rank(comm,&rank);

1340:   dd->dim         = 2;
1341:   PetscMalloc(dof*sizeof(char*),&dd->fieldname);
1342:   PetscMemzero(dd->fieldname,dof*sizeof(char*));

1344:   if (m != PETSC_DECIDE) {
1345:     if (m < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in X direction: %D",m);
1346:     else if (m > size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in X direction: %D %d",m,size);
1347:   }
1348:   if (n != PETSC_DECIDE) {
1349:     if (n < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Y direction: %D",n);
1350:     else if (n > size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Y direction: %D %d",n,size);
1351:   }

1353:   if (m == PETSC_DECIDE || n == PETSC_DECIDE) {
1354:     if (n != PETSC_DECIDE) {
1355:       m = size/n;
1356:     } else if (m != PETSC_DECIDE) {
1357:       n = size/m;
1358:     } else {
1359:       /* try for squarish distribution */
1360:       m = (PetscInt)(0.5 + sqrt(((double)M)*((double)size)/((double)N)));
1361:       if (!m) m = 1;
1362:       while (m > 0) {
1363:         n = size/m;
1364:         if (m*n == size) break;
1365:         m--;
1366:       }
1367:       if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;}
1368:     }
1369:     if (m*n != size) SETERRQ(comm,PETSC_ERR_PLIB,"Unable to create partition, check the size of the communicator and input m and n ");
1370:   } else if (m*n != size) SETERRQ(comm,PETSC_ERR_ARG_OUTOFRANGE,"Given Bad partition");

1372:   if (M < m) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Partition in x direction is too fine! %D %D",M,m);
1373:   if (N < n) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Partition in y direction is too fine! %D %D",N,n);

1375:   /* 
1376:      Determine locally owned region 
1377:      xs is the first local node number, x is the number of local nodes 
1378:   */
1379:   if (!lx) {
1380:     PetscMalloc(m*sizeof(PetscInt), &dd->lx);
1381:     lx = dd->lx;
1382:     for (i=0; i<m; i++) {
1383:       lx[i] = M/m + ((M % m) > i);
1384:     }
1385:   }
1386:   x  = lx[rank % m];
1387:   xs = 0;
1388:   for (i=0; i<(rank % m); i++) {
1389:     xs += lx[i];
1390:   }
1391: #if defined(PETSC_USE_DEBUG)
1392:   left = xs;
1393:   for (i=(rank % m); i<m; i++) {
1394:     left += lx[i];
1395:   }
1396:   if (left != M) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Sum of lx across processors not equal to M: %D %D",left,M);
1397: #endif

1399:   /* 
1400:      Determine locally owned region 
1401:      ys is the first local node number, y is the number of local nodes 
1402:   */
1403:   if (!ly) {
1404:     PetscMalloc(n*sizeof(PetscInt), &dd->ly);
1405:     ly = dd->ly;
1406:     for (i=0; i<n; i++) {
1407:       ly[i] = N/n + ((N % n) > i);
1408:     }
1409:   }
1410:   y  = ly[rank/m];
1411:   ys = 0;
1412:   for (i=0; i<(rank/m); i++) {
1413:     ys += ly[i];
1414:   }
1415: #if defined(PETSC_USE_DEBUG)
1416:   left = ys;
1417:   for (i=(rank/m); i<n; i++) {
1418:     left += ly[i];
1419:   }
1420:   if (left != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Sum of ly across processors not equal to N: %D %D",left,N);
1421: #endif

1423:   /*
1424:    check if the scatter requires more than one process neighbor or wraps around
1425:    the domain more than once
1426:   */
1427:   if ((x < s) && ((m > 1) || (bx == DMDA_BOUNDARY_PERIODIC))) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local x-width of domain x %D is smaller than stencil width s %D",x,s);
1428:   if ((y < s) && ((n > 1) || (by == DMDA_BOUNDARY_PERIODIC))) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local y-width of domain y %D is smaller than stencil width s %D",y,s);
1429:   xe = xs + x;
1430:   ye = ys + y;

1432:   /* determine ghost region (Xs) and region scattered into (IXs)  */
1433:   /* Assume No Periodicity */
1434:   if (xs-s > 0) { Xs = xs - s; IXs = xs - s; } else { Xs = 0; IXs = 0; }
1435:   if (xe+s <= M) { Xe = xe + s; IXe = xe + s; } else { Xe = M; IXe = M; }
1436:   if (ys-s > 0) { Ys = ys - s; IYs = ys - s; } else { Ys = 0; IYs = 0; }
1437:   if (ye+s <= N) { Ye = ye + s; IYe = ye + s; } else { Ye = N; IYe = N; }

1439:   /* fix for periodicity/ghosted */
1440:   if (bx) { Xs = xs - s; Xe = xe + s; }
1441:   if (bx == DMDA_BOUNDARY_PERIODIC) { IXs = xs - s; IXe = xe + s; }
1442:   if (by) { Ys = ys - s; Ye = ye + s; }
1443:   if (by == DMDA_BOUNDARY_PERIODIC) { IYs = ys - s; IYe = ye + s; }

1445:   /* Resize all X parameters to reflect w */
1446:   s_x = s;
1447:   s_y = s;

1449:   /* determine starting point of each processor */
1450:   nn    = x*y;
1451:   PetscMalloc2(size+1,PetscInt,&bases,size,PetscInt,&ldims);
1452:   MPI_Allgather(&nn,1,MPIU_INT,ldims,1,MPIU_INT,comm);
1453:   bases[0] = 0;
1454:   for (i=1; i<=size; i++) {
1455:     bases[i] = ldims[i-1];
1456:   }
1457:   for (i=1; i<=size; i++) {
1458:     bases[i] += bases[i-1];
1459:   }
1460:   base = bases[rank]*dof;

1462:   /* allocate the base parallel and sequential vectors */
1463:   dd->Nlocal = x*y*dof;
1464:   VecCreateMPIWithArray(comm,dof,dd->Nlocal,PETSC_DECIDE,0,&global);
1465:   dd->nlocal = (Xe-Xs)*(Ye-Ys)*dof;
1466:   VecCreateSeqWithArray(PETSC_COMM_SELF,dof,dd->nlocal,0,&local);

1468:   /* generate appropriate vector scatters */
1469:   /* local to global inserts non-ghost point region into global */
1470:   VecGetOwnershipRange(global,&start,&end);
1471:   ISCreateStride(comm,x*y*dof,start,1,&to);

1473:   count = x*y;
1474:   PetscMalloc(x*y*sizeof(PetscInt),&idx);
1475:   left = xs - Xs; right = left + x;
1476:   down = ys - Ys; up = down + y;
1477:   count = 0;
1478:   for (i=down; i<up; i++) {
1479:     for (j=left; j<right; j++) {
1480:       idx[count++] = i*(Xe-Xs) + j;
1481:     }
1482:   }

1484:   ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&from);
1485:   VecScatterCreate(local,from,global,to,&ltog);
1486:   PetscLogObjectParent(dd,ltog);
1487:   ISDestroy(&from);
1488:   ISDestroy(&to);

1490:   /* global to local must include ghost points within the domain,
1491:      but not ghost points outside the domain that aren't periodic */
1492:   if (stencil_type == DMDA_STENCIL_BOX) {
1493:     count = (IXe-IXs)*(IYe-IYs);
1494:     PetscMalloc(count*sizeof(PetscInt),&idx);

1496:     left = IXs - Xs; right = left + (IXe-IXs);
1497:     down = IYs - Ys; up = down + (IYe-IYs);
1498:     count = 0;
1499:     for (i=down; i<up; i++) {
1500:       for (j=left; j<right; j++) {
1501:         idx[count++] = j + i*(Xe-Xs);
1502:       }
1503:     }
1504:     ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);

1506:   } else {
1507:     /* must drop into cross shape region */
1508:     /*       ---------|
1509:             |  top    |
1510:          |---         ---| up
1511:          |   middle      |
1512:          |               |
1513:          ----         ---- down
1514:             | bottom  |
1515:             -----------
1516:          Xs xs        xe Xe */
1517:     count = (ys-IYs)*x + y*(IXe-IXs) + (IYe-ye)*x;
1518:     PetscMalloc(count*sizeof(PetscInt),&idx);

1520:     left = xs - Xs; right = left + x;
1521:     down = ys - Ys; up = down + y;
1522:     count = 0;
1523:     /* bottom */
1524:     for (i=(IYs-Ys); i<down; i++) {
1525:       for (j=left; j<right; j++) {
1526:         idx[count++] = j + i*(Xe-Xs);
1527:       }
1528:     }
1529:     /* middle */
1530:     for (i=down; i<up; i++) {
1531:       for (j=(IXs-Xs); j<(IXe-Xs); j++) {
1532:         idx[count++] = j + i*(Xe-Xs);
1533:       }
1534:     }
1535:     /* top */
1536:     for (i=up; i<up+IYe-ye; i++) {
1537:       for (j=left; j<right; j++) {
1538:         idx[count++] = j + i*(Xe-Xs);
1539:       }
1540:     }
1541:     ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);
1542:   }


1545:   /* determine who lies on each side of us stored in    n6 n7 n8
1546:                                                         n3    n5
1547:                                                         n0 n1 n2
1548:   */

1550:   /* Assume the Non-Periodic Case */
1551:   n1 = rank - m;
1552:   if (rank % m) {
1553:     n0 = n1 - 1;
1554:   } else {
1555:     n0 = -1;
1556:   }
1557:   if ((rank+1) % m) {
1558:     n2 = n1 + 1;
1559:     n5 = rank + 1;
1560:     n8 = rank + m + 1; if (n8 >= m*n) n8 = -1;
1561:   } else {
1562:     n2 = -1; n5 = -1; n8 = -1;
1563:   }
1564:   if (rank % m) {
1565:     n3 = rank - 1;
1566:     n6 = n3 + m; if (n6 >= m*n) n6 = -1;
1567:   } else {
1568:     n3 = -1; n6 = -1;
1569:   }
1570:   n7 = rank + m; if (n7 >= m*n) n7 = -1;

1572:   if (bx == DMDA_BOUNDARY_PERIODIC && by == DMDA_BOUNDARY_PERIODIC) {
1573:   /* Modify for Periodic Cases */
1574:     /* Handle all four corners */
1575:     if ((n6 < 0) && (n7 < 0) && (n3 < 0)) n6 = m-1;
1576:     if ((n8 < 0) && (n7 < 0) && (n5 < 0)) n8 = 0;
1577:     if ((n2 < 0) && (n5 < 0) && (n1 < 0)) n2 = size-m;
1578:     if ((n0 < 0) && (n3 < 0) && (n1 < 0)) n0 = size-1;

1580:     /* Handle Top and Bottom Sides */
1581:     if (n1 < 0) n1 = rank + m * (n-1);
1582:     if (n7 < 0) n7 = rank - m * (n-1);
1583:     if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
1584:     if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1;
1585:     if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
1586:     if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1;

1588:     /* Handle Left and Right Sides */
1589:     if (n3 < 0) n3 = rank + (m-1);
1590:     if (n5 < 0) n5 = rank - (m-1);
1591:     if ((n1 >= 0) && (n0 < 0)) n0 = rank-1;
1592:     if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1;
1593:     if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1;
1594:     if ((n7 >= 0) && (n8 < 0)) n8 = rank+1;
1595:   } else if (by == DMDA_BOUNDARY_PERIODIC) {  /* Handle Top and Bottom Sides */
1596:     if (n1 < 0) n1 = rank + m * (n-1);
1597:     if (n7 < 0) n7 = rank - m * (n-1);
1598:     if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
1599:     if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1;
1600:     if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
1601:     if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1;
1602:   } else if (bx == DMDA_BOUNDARY_PERIODIC) { /* Handle Left and Right Sides */
1603:     if (n3 < 0) n3 = rank + (m-1);
1604:     if (n5 < 0) n5 = rank - (m-1);
1605:     if ((n1 >= 0) && (n0 < 0)) n0 = rank-1;
1606:     if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1;
1607:     if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1;
1608:     if ((n7 >= 0) && (n8 < 0)) n8 = rank+1;
1609:   }

1611:   PetscMalloc(9*sizeof(PetscInt),&dd->neighbors);
1612:   dd->neighbors[0] = n0;
1613:   dd->neighbors[1] = n1;
1614:   dd->neighbors[2] = n2;
1615:   dd->neighbors[3] = n3;
1616:   dd->neighbors[4] = rank;
1617:   dd->neighbors[5] = n5;
1618:   dd->neighbors[6] = n6;
1619:   dd->neighbors[7] = n7;
1620:   dd->neighbors[8] = n8;

1622:   if (stencil_type == DMDA_STENCIL_STAR) {
1623:     /* save corner processor numbers */
1624:     sn0 = n0; sn2 = n2; sn6 = n6; sn8 = n8;
1625:     n0 = n2 = n6 = n8 = -1;
1626:   }

1628:   PetscMalloc((Xe-Xs)*(Ye-Ys)*sizeof(PetscInt),&idx);
1629:   PetscLogObjectMemory(da,(Xe-Xs)*(Ye-Ys)*sizeof(PetscInt));

1631:   nn = 0;
1632:   xbase = bases[rank];
1633:   for (i=1; i<=s_y; i++) {
1634:     if (n0 >= 0) { /* left below */
1635:       x_t = lx[n0 % m];
1636:       y_t = ly[(n0/m)];
1637:       s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x;
1638:       for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1639:     }
1640:     if (n1 >= 0) { /* directly below */
1641:       x_t = x;
1642:       y_t = ly[(n1/m)];
1643:       s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t;
1644:       for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1645:     }
1646:     if (n2 >= 0) { /* right below */
1647:       x_t = lx[n2 % m];
1648:       y_t = ly[(n2/m)];
1649:       s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t;
1650:       for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1651:     }
1652:   }

1654:   for (i=0; i<y; i++) {
1655:     if (n3 >= 0) { /* directly left */
1656:       x_t = lx[n3 % m];
1657:       /* y_t = y; */
1658:       s_t = bases[n3] + (i+1)*x_t - s_x;
1659:       for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1660:     }

1662:     for (j=0; j<x; j++) { idx[nn++] = xbase++; } /* interior */

1664:     if (n5 >= 0) { /* directly right */
1665:       x_t = lx[n5 % m];
1666:       /* y_t = y; */
1667:       s_t = bases[n5] + (i)*x_t;
1668:       for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1669:     }
1670:   }

1672:   for (i=1; i<=s_y; i++) {
1673:     if (n6 >= 0) { /* left above */
1674:       x_t = lx[n6 % m];
1675:       /* y_t = ly[(n6/m)]; */
1676:       s_t = bases[n6] + (i)*x_t - s_x;
1677:       for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1678:     }
1679:     if (n7 >= 0) { /* directly above */
1680:       x_t = x;
1681:       /* y_t = ly[(n7/m)]; */
1682:       s_t = bases[n7] + (i-1)*x_t;
1683:       for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1684:     }
1685:     if (n8 >= 0) { /* right above */
1686:       x_t = lx[n8 % m];
1687:       /* y_t = ly[(n8/m)]; */
1688:       s_t = bases[n8] + (i-1)*x_t;
1689:       for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1690:     }
1691:   }

1693:   ISCreateBlock(comm,dof,nn,idx,PETSC_COPY_VALUES,&from);
1694:   VecScatterCreate(global,from,local,to,&gtol);
1695:   PetscLogObjectParent(da,gtol);
1696:   ISDestroy(&to);
1697:   ISDestroy(&from);

1699:   if (stencil_type == DMDA_STENCIL_STAR) {
1700:     n0 = sn0; n2 = sn2; n6 = sn6; n8 = sn8;
1701:   }

1703:   if ((stencil_type == DMDA_STENCIL_STAR) ||
1704:       (bx && bx != DMDA_BOUNDARY_PERIODIC) ||
1705:       (by && by != DMDA_BOUNDARY_PERIODIC)) {
1706:     /*
1707:         Recompute the local to global mappings, this time keeping the 
1708:       information about the cross corner processor numbers and any ghosted
1709:       but not periodic indices.
1710:     */
1711:     nn = 0;
1712:     xbase = bases[rank];
1713:     for (i=1; i<=s_y; i++) {
1714:       if (n0 >= 0) { /* left below */
1715:         x_t = lx[n0 % m];
1716:         y_t = ly[(n0/m)];
1717:         s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x;
1718:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1719:       } else if (xs-Xs > 0 && ys-Ys > 0) {
1720:         for (j=0; j<s_x; j++) { idx[nn++] = -1;}
1721:       }
1722:       if (n1 >= 0) { /* directly below */
1723:         x_t = x;
1724:         y_t = ly[(n1/m)];
1725:         s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t;
1726:         for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1727:       } else if (ys-Ys > 0) {
1728:         for (j=0; j<x; j++) { idx[nn++] = -1;}
1729:       }
1730:       if (n2 >= 0) { /* right below */
1731:         x_t = lx[n2 % m];
1732:         y_t = ly[(n2/m)];
1733:         s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t;
1734:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1735:       } else if (Xe-xe> 0 && ys-Ys > 0) {
1736:         for (j=0; j<s_x; j++) { idx[nn++] = -1;}
1737:       }
1738:     }

1740:     for (i=0; i<y; i++) {
1741:       if (n3 >= 0) { /* directly left */
1742:         x_t = lx[n3 % m];
1743:         /* y_t = y; */
1744:         s_t = bases[n3] + (i+1)*x_t - s_x;
1745:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1746:       } else if (xs-Xs > 0) {
1747:         for (j=0; j<s_x; j++) { idx[nn++] = -1;}
1748:       }

1750:       for (j=0; j<x; j++) { idx[nn++] = xbase++; } /* interior */

1752:       if (n5 >= 0) { /* directly right */
1753:         x_t = lx[n5 % m];
1754:         /* y_t = y; */
1755:         s_t = bases[n5] + (i)*x_t;
1756:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1757:       } else if (Xe-xe > 0) {
1758:         for (j=0; j<s_x; j++) { idx[nn++] = -1;}
1759:       }
1760:     }

1762:     for (i=1; i<=s_y; i++) {
1763:       if (n6 >= 0) { /* left above */
1764:         x_t = lx[n6 % m];
1765:         /* y_t = ly[(n6/m)]; */
1766:         s_t = bases[n6] + (i)*x_t - s_x;
1767:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1768:       } else if (xs-Xs > 0 && Ye-ye > 0) {
1769:         for (j=0; j<s_x; j++) { idx[nn++] = -1;}
1770:       }
1771:       if (n7 >= 0) { /* directly above */
1772:         x_t = x;
1773:         /* y_t = ly[(n7/m)]; */
1774:         s_t = bases[n7] + (i-1)*x_t;
1775:         for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1776:       } else if (Ye-ye > 0) {
1777:         for (j=0; j<x; j++) { idx[nn++] = -1;}
1778:       }
1779:       if (n8 >= 0) { /* right above */
1780:         x_t = lx[n8 % m];
1781:         /* y_t = ly[(n8/m)]; */
1782:         s_t = bases[n8] + (i-1)*x_t;
1783:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1784:       } else if (Xe-xe > 0 && Ye-ye > 0) {
1785:         for (j=0; j<s_x; j++) { idx[nn++] = -1;}
1786:       }
1787:     }
1788:   }
1789:   /*
1790:      Set the local to global ordering in the global vector, this allows use
1791:      of VecSetValuesLocal().
1792:   */
1793:   ISCreateBlock(comm,dof,nn,idx,PETSC_OWN_POINTER,&ltogis);
1794:   PetscMalloc(nn*dof*sizeof(PetscInt),&idx_cpy);
1795:   PetscLogObjectMemory(da,nn*dof*sizeof(PetscInt));
1796:   ISGetIndices(ltogis, &idx_full);
1797:   PetscMemcpy(idx_cpy,idx_full,nn*dof*sizeof(PetscInt));
1798:   ISRestoreIndices(ltogis, &idx_full);
1799:   ISLocalToGlobalMappingCreateIS(ltogis,&da->ltogmap);
1800:   PetscLogObjectParent(da,da->ltogmap);
1801:   ISDestroy(&ltogis);
1802:   ISLocalToGlobalMappingBlock(da->ltogmap,dd->w,&da->ltogmapb);
1803:   PetscLogObjectParent(da,da->ltogmap);

1805:   PetscFree2(bases,ldims);
1806:   dd->m  = m;  dd->n  = n;
1807:   /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */
1808:   dd->xs = xs*dof; dd->xe = xe*dof; dd->ys = ys; dd->ye = ye; dd->zs = 0; dd->ze = 1;
1809:   dd->Xs = Xs*dof; dd->Xe = Xe*dof; dd->Ys = Ys; dd->Ye = Ye; dd->Zs = 0; dd->Ze = 1;

1811:   VecDestroy(&local);
1812:   VecDestroy(&global);

1814:   dd->gtol      = gtol;
1815:   dd->ltog      = ltog;
1816:   dd->idx       = idx_cpy;
1817:   dd->Nl        = nn*dof;
1818:   dd->base      = base;
1819:   da->ops->view = DMView_DA_2d;
1820:   dd->ltol = PETSC_NULL;
1821:   dd->ao   = PETSC_NULL;

1823:   return(0);
1824: }

1828: /*@C
1829:    DMDACreate2d -  Creates an object that will manage the communication of  two-dimensional 
1830:    regular array data that is distributed across some processors.

1832:    Collective on MPI_Comm

1834:    Input Parameters:
1835: +  comm - MPI communicator
1836: .  bx,by - type of ghost nodes the array have. 
1837:          Use one of DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_PERIODIC.
1838: .  stencil_type - stencil type.  Use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.
1839: .  M,N - global dimension in each direction of the array (use -M and or -N to indicate that it may be set to a different value 
1840:             from the command line with -da_grid_x <M> -da_grid_y <N>)
1841: .  m,n - corresponding number of processors in each dimension 
1842:          (or PETSC_DECIDE to have calculated)
1843: .  dof - number of degrees of freedom per node
1844: .  s - stencil width
1845: -  lx, ly - arrays containing the number of nodes in each cell along
1846:            the x and y coordinates, or PETSC_NULL. If non-null, these
1847:            must be of length as m and n, and the corresponding
1848:            m and n cannot be PETSC_DECIDE. The sum of the lx[] entries
1849:            must be M, and the sum of the ly[] entries must be N.

1851:    Output Parameter:
1852: .  da - the resulting distributed array object

1854:    Options Database Key:
1855: +  -da_view - Calls DMView() at the conclusion of DMDACreate2d()
1856: .  -da_grid_x <nx> - number of grid points in x direction, if M < 0
1857: .  -da_grid_y <ny> - number of grid points in y direction, if N < 0
1858: .  -da_processors_x <nx> - number of processors in x direction
1859: .  -da_processors_y <ny> - number of processors in y direction
1860: .  -da_refine_x <rx> - refinement ratio in x direction
1861: .  -da_refine_y <ry> - refinement ratio in y direction
1862: -  -da_refine <n> - refine the DMDA n times before creating, if M or N < 0


1865:    Level: beginner

1867:    Notes:
1868:    The stencil type DMDA_STENCIL_STAR with width 1 corresponds to the 
1869:    standard 5-pt stencil, while DMDA_STENCIL_BOX with width 1 denotes
1870:    the standard 9-pt stencil.

1872:    The array data itself is NOT stored in the DMDA, it is stored in Vec objects;
1873:    The appropriate vector objects can be obtained with calls to DMCreateGlobalVector()
1874:    and DMCreateLocalVector() and calls to VecDuplicate() if more are needed.

1876: .keywords: distributed array, create, two-dimensional

1878: .seealso: DMDestroy(), DMView(), DMDACreate1d(), DMDACreate3d(), DMGlobalToLocalBegin(), DMDAGetRefinementFactor(),
1879:           DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMDALocalToLocalBegin(), DMDALocalToLocalEnd(), DMDASetRefinementFactor(),
1880:           DMDAGetInfo(), DMCreateGlobalVector(), DMCreateLocalVector(), DMDACreateNaturalVector(), DMLoad(), DMDAGetOwnershipRanges()

1882: @*/

1884: PetscErrorCode  DMDACreate2d(MPI_Comm comm,DMDABoundaryType bx,DMDABoundaryType by,DMDAStencilType stencil_type,
1885:                           PetscInt M,PetscInt N,PetscInt m,PetscInt n,PetscInt dof,PetscInt s,const PetscInt lx[],const PetscInt ly[],DM *da)
1886: {

1890:   DMDACreate(comm, da);
1891:   DMDASetDim(*da, 2);
1892:   DMDASetSizes(*da, M, N, 1);
1893:   DMDASetNumProcs(*da, m, n, PETSC_DECIDE);
1894:   DMDASetBoundaryType(*da, bx, by, DMDA_BOUNDARY_NONE);
1895:   DMDASetDof(*da, dof);
1896:   DMDASetStencilType(*da, stencil_type);
1897:   DMDASetStencilWidth(*da, s);
1898:   DMDASetOwnershipRanges(*da, lx, ly, PETSC_NULL);
1899:   /* This violates the behavior for other classes, but right now users expect negative dimensions to be handled this way */
1900:   DMSetFromOptions(*da);
1901:   DMSetUp(*da);
1902:   DMView_DA_Private(*da);
1903:   return(0);
1904: }