自学教程：C++ DMGlobalToLocalBegin函数代码示例

static PetscErrorCode FormFunction_All(SNES snes,Vec X,Vec F,void *ctx){  User           user = (User)ctx;  DM             dau,dak;  DMDALocalInfo  infou,infok;  PetscScalar    *u,*k;  PetscScalar    *fu,*fk;  PetscErrorCode ierr;  Vec            Uloc,Kloc,Fu,Fk;  PetscFunctionBeginUser;  ierr = DMCompositeGetEntries(user->pack,&dau,&dak);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(dau,&infou);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(dak,&infok);CHKERRQ(ierr);  ierr = DMCompositeGetLocalVectors(user->pack,&Uloc,&Kloc);CHKERRQ(ierr);  switch (user->ptype) {  case 0:    ierr = DMGlobalToLocalBegin(dau,X,INSERT_VALUES,Uloc);CHKERRQ(ierr);    ierr = DMGlobalToLocalEnd  (dau,X,INSERT_VALUES,Uloc);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dau,Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dak,user->Kloc,&k);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dau,F,&fu);CHKERRQ(ierr);    ierr = FormFunctionLocal_U(user,&infou,u,k,fu);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dau,F,&fu);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dau,Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dak,user->Kloc,&k);CHKERRQ(ierr);    break;  case 1:    ierr = DMGlobalToLocalBegin(dak,X,INSERT_VALUES,Kloc);CHKERRQ(ierr);    ierr = DMGlobalToLocalEnd  (dak,X,INSERT_VALUES,Kloc);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dau,user->Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dak,Kloc,&k);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dak,F,&fk);CHKERRQ(ierr);    ierr = FormFunctionLocal_K(user,&infok,u,k,fk);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dak,F,&fk);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dau,user->Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dak,Kloc,&k);CHKERRQ(ierr);    break;  case 2:    ierr = DMCompositeScatter(user->pack,X,Uloc,Kloc);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dau,Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dak,Kloc,&k);CHKERRQ(ierr);    ierr = DMCompositeGetAccess(user->pack,F,&Fu,&Fk);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dau,Fu,&fu);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dak,Fk,&fk);CHKERRQ(ierr);    ierr = FormFunctionLocal_U(user,&infou,u,k,fu);CHKERRQ(ierr);    ierr = FormFunctionLocal_K(user,&infok,u,k,fk);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dau,Fu,&fu);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dak,Fk,&fk);CHKERRQ(ierr);    ierr = DMCompositeRestoreAccess(user->pack,F,&Fu,&Fk);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dau,Uloc,&u);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dak,Kloc,&k);CHKERRQ(ierr);    break;  }  ierr = DMCompositeRestoreLocalVectors(user->pack,&Uloc,&Kloc);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx){  PetscErrorCode ierr;  AppCtx         *user=(AppCtx*)ctx;  DM             cda;  DMDACoor2d     **coors;  PetscScalar    **p,**f,**pdot;  PetscInt       i,j;  PetscInt       xs,ys,xm,ym,M,N;  Vec            localX,gc,localXdot;  PetscScalar    p_adv1,p_adv2,p_diff;  PetscFunctionBeginUser;  ierr = DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);  ierr = DMGetCoordinateDM(user->da,&cda);CHKERRQ(ierr);  ierr = DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);  ierr = DMGetLocalVector(user->da,&localX);CHKERRQ(ierr);  ierr = DMGetLocalVector(user->da,&localXdot);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(user->da,Xdot,INSERT_VALUES,localXdot);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(user->da,Xdot,INSERT_VALUES,localXdot);CHKERRQ(ierr);  ierr = DMGetCoordinatesLocal(user->da,&gc);CHKERRQ(ierr);  ierr = DMDAVecGetArrayRead(cda,gc,&coors);CHKERRQ(ierr);  ierr = DMDAVecGetArrayRead(user->da,localX,&p);CHKERRQ(ierr);  ierr = DMDAVecGetArrayRead(user->da,localXdot,&pdot);CHKERRQ(ierr);  ierr = DMDAVecGetArray(user->da,F,&f);CHKERRQ(ierr);  PetscScalar diffuse1,gamma;  gamma = user->D*user->ws/(2*user->H);  diffuse1 = user->lambda*user->lambda*user->q/(user->lambda*gamma+1)*(1.0 - PetscExpScalar(-t*(gamma+1.0)/user->lambda));  user->disper_coe = user->ws*user->ws/(4*user->H*user->H)*diffuse1;  for (i=xs; i < xs+xm; i++) {    for (j=ys; j < ys+ym; j++) {      ierr = adv1(p,coors[j][i].y,i,j,M,&p_adv1,user);CHKERRQ(ierr);      ierr = adv2(p,coors[j][i].x,coors[j][i].y,i,j,N,&p_adv2,user);CHKERRQ(ierr);      ierr = diffuse(p,i,j,t,&p_diff,user);CHKERRQ(ierr);      f[j][i] = -p_adv1 - p_adv2  + p_diff - pdot[j][i];    }  }  ierr = DMDAVecRestoreArrayRead(user->da,localX,&p);CHKERRQ(ierr);  ierr = DMDAVecRestoreArrayRead(user->da,localX,&pdot);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(user->da,&localX);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(user->da,&localXdot);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(user->da,F,&f);CHKERRQ(ierr);  ierr = DMDAVecRestoreArrayRead(cda,gc,&coors);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode FormResidual(TS ts,PetscReal t,Vec U,Vec Udot,Vec R,void *user){  PetscFunctionBegin;  PetscErrorCode ierr;  DMDALocalInfo    info;  DM               dm;  Vec              localU,localUdot,localR; // local versions  Field          **h,**hdot,**r;  /* get the da from the snes */  ierr = TSGetDM(ts, &dm);CHKERRQ(ierr);  /* handle the vec U */  ierr = DMGetLocalVector(dm,&localU);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(dm,U,INSERT_VALUES,localU);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(dm,U,INSERT_VALUES,localU);CHKERRQ(ierr);  /* handle the vec Udot */  ierr = DMGetLocalVector(dm,&localUdot);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(dm,Udot,INSERT_VALUES,localUdot);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(dm,Udot,INSERT_VALUES,localUdot);CHKERRQ(ierr);  /* handle the vec R */  ierr = DMGetLocalVector(dm,&localR);CHKERRQ(ierr);  ierr = VecZeroEntries(localR);CHKERRQ(ierr);  /* Get the arrays from the vectors */  ierr = DMIGAVecGetArray(dm,localU,&h);CHKERRQ(ierr);  ierr = DMIGAVecGetArray(dm,localUdot,&hdot);CHKERRQ(ierr);  ierr = DMIGAVecGetArray(dm,localR,&r);CHKERRQ(ierr);  /* Grab the local info and call the local residual routine */  ierr = DMIGAGetLocalInfo(dm,&info);CHKERRQ(ierr);CHKERRQ(ierr);  ierr = FormResidualLocal(&info,t,h,hdot,r,(AppCtx *) user);CHKERRQ(ierr);  /* Restore the arrays */  ierr = DMIGAVecRestoreArray(dm,localR,&r);CHKERRQ(ierr);  ierr = DMIGAVecRestoreArray(dm,localUdot,&hdot);CHKERRQ(ierr);  ierr = DMIGAVecRestoreArray(dm,localU,&h);CHKERRQ(ierr);  /* Add contributions back to global R */  ierr = VecZeroEntries(R);CHKERRQ(ierr);  ierr = DMLocalToGlobalBegin(dm,localR,ADD_VALUES,R);CHKERRQ(ierr); // this one adds the values  ierr = DMLocalToGlobalEnd(dm,localR,ADD_VALUES,R);CHKERRQ(ierr);  /* Restore the local vectors */  ierr = DMRestoreLocalVector(dm,&localU);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm,&localUdot);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm,&localR);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode FormTangent(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flag,void *ctx){  PetscFunctionBegin;  PetscErrorCode ierr;  SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "FormTangent not implemented, use -snes_mf");  DMDALocalInfo    info;  DM               da_dof;  Vec              localU,localUdot; // local versions  Field          **h,**hdot;  /* get the da from the snes */  ierr = TSGetDM(ts,(DM*)&da_dof);CHKERRQ(ierr);  /* handle the vec U */  ierr = DMGetLocalVector(da_dof,&localU);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(da_dof,U,INSERT_VALUES,localU);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da_dof,U,INSERT_VALUES,localU);CHKERRQ(ierr);  /* handle the vec Udot */  ierr = DMGetLocalVector(da_dof,&localUdot);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(da_dof,Udot,INSERT_VALUES,localUdot);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da_dof,Udot,INSERT_VALUES,localUdot);CHKERRQ(ierr);  /* Get the arrays from the vectors */  ierr = DMDAVecGetArray(da_dof,localU,&h);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da_dof,localUdot,&hdot);CHKERRQ(ierr);  /* Grab the local info and call the local tangent routine */  ierr = DMDAGetLocalInfo(da_dof,&info);CHKERRQ(ierr);CHKERRQ(ierr);  ierr = MatZeroEntries(*B);CHKERRQ(ierr); // pre-zero the matrix  ierr = FormTangentLocal(&info,t,h,hdot,shift,B,(AppCtx *) ctx);CHKERRQ(ierr);  ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);  ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);  if (*A != *B) { // then we could be matrix free    ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);    ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);  }  *flag = SAME_NONZERO_PATTERN; /* the sparsity pattern does not change */  ierr = DMDAVecRestoreArray(da_dof,localUdot,&hdot);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da_dof,localU,&h);CHKERRQ(ierr);  /* Restore the arrays and local vectors */  ierr = DMRestoreLocalVector(da_dof,&localU);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da_dof,&localUdot);CHKERRQ(ierr);  PetscFunctionReturn(0);}

/*@C  DMPlexProjectField - This projects the given function of the fields into the function space provided.  Input Parameters:+ dm      - The DM. U       - The input field vector. funcs   - The functions to evaluate, one per field- mode    - The insertion mode for values  Output Parameter:. X       - The output vector  Level: developer.seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff()@*/PetscErrorCode DMPlexProjectField(DM dm, Vec U, void (**funcs)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal [], PetscScalar []), InsertMode mode, Vec X){  Vec            localX, localU;  PetscErrorCode ierr;  PetscFunctionBegin;  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);  ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm, &localU);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(dm, U, INSERT_VALUES, localU);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(dm, U, INSERT_VALUES, localU);CHKERRQ(ierr);  ierr = DMPlexProjectFieldLocal(dm, localU, funcs, mode, localX);CHKERRQ(ierr);  ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr);  if (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES) {    Mat cMat;    ierr = DMGetDefaultConstraints(dm, NULL, &cMat);CHKERRQ(ierr);    if (cMat) {      ierr = DMGlobalToLocalSolve(dm, localX, X);CHKERRQ(ierr);    }  }  ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm, &localU);CHKERRQ(ierr);  PetscFunctionReturn(0);}

/**   Form the time step function for the mass term   @param Vec in:       (y^(n+1))   out @param Vec out:  My^(n+1) - M y^(n)    if timestep = 1, backward euler                                              will do the BDF2**/PetscErrorCode FormMassTimeStepFunction(User user, Algebra algebra, Vec in, Vec out, PetscBool rebuild){  PetscErrorCode  ierr;  PetscMPIInt     rank;  Vec            inLocal;  PetscFunctionBegin;  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);  ierr = VecSet(out, 0.0);CHKERRQ(ierr);  ierr = DMGetLocalVector(user->dm, &inLocal);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(user->dm, in, INSERT_VALUES, inLocal);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(user->dm, in, INSERT_VALUES, inLocal);CHKERRQ(ierr);  ierr = ApplyBC(user->dm, user->current_time, inLocal, user);CHKERRQ(ierr);  ierr = CaculateLocalMassFunction(user->dm, inLocal, out, user);CHKERRQ(ierr);  ierr = DMLocalToGlobalBegin(user->dm, inLocal, INSERT_VALUES, in);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(user->dm, inLocal, INSERT_VALUES, in);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(user->dm, &inLocal);CHKERRQ(ierr);  PetscFunctionReturn(0);}

/*  This function should eventually replace:    DMDAComputeFunction() and DMDAComputeFunction1() */static PetscErrorCode SNESComputeObjective_DMDA(SNES snes,Vec X,PetscReal *ob,void *ctx){  PetscErrorCode ierr;  DM             dm;  DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;  DMDALocalInfo  info;  Vec            Xloc;  void           *x;  PetscFunctionBegin;  PetscValidHeaderSpecific(snes,SNES_CLASSID,1);  PetscValidHeaderSpecific(X,VEC_CLASSID,2);  PetscValidPointer(ob,3);  if (!dmdasnes->objectivelocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");  ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);  ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);  CHKMEMQ;  ierr = (*dmdasnes->objectivelocal)(&info,x,ob,dmdasnes->objectivelocalctx);CHKERRQ(ierr);  CHKMEMQ;  ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);  PetscFunctionReturn(0);}

static PetscErrorCode SNESComputeFunction_DMLocal(SNES snes,Vec X,Vec F,void *ctx){  PetscErrorCode ierr;  DM             dm;  DMSNES_Local   *dmlocalsnes = (DMSNES_Local*)ctx;  Vec            Xloc,Floc;  PetscFunctionBegin;  PetscValidHeaderSpecific(snes,SNES_CLASSID,1);  PetscValidHeaderSpecific(X,VEC_CLASSID,2);  PetscValidHeaderSpecific(F,VEC_CLASSID,3);  ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm,&Floc);CHKERRQ(ierr);  ierr = VecZeroEntries(Xloc);CHKERRQ(ierr);  ierr = VecZeroEntries(Floc);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  CHKMEMQ;  ierr = (*dmlocalsnes->residuallocal)(dm,Xloc,Floc,dmlocalsnes->residuallocalctx);CHKERRQ(ierr);  CHKMEMQ;  ierr = VecZeroEntries(F);CHKERRQ(ierr);  ierr = DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm,&Floc);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);  PetscFunctionReturn(0);}

/*   Applies the second order centered difference diffusion operator on a one dimensional periodic domain*/static PetscErrorCode FormDiffusionFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr){  User              user = (User)ptr;  PetscErrorCode    ierr;  PetscScalar       **f;  const PetscScalar **x;  DM                dm;  PetscInt          i,xs,xm,j,dof;  Vec               Xlocal;  PetscReal         idx;  PetscFunctionBeginUser;  ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);  ierr = DMDAGetInfo(dm,NULL,NULL,NULL,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm,&Xlocal);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xlocal);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xlocal);CHKERRQ(ierr);  ierr = DMDAVecGetArrayDOFRead(dm,Xlocal,&x);CHKERRQ(ierr);  ierr = DMDAVecGetArrayDOF(dm,F,&f);CHKERRQ(ierr);  ierr = DMDAGetCorners(dm,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);  idx = 1.0*user->diffus/user->dx;  for (i=xs; i<xs+xm; i++) {    for (j=0; j<dof; j++) {      f[i][j] += idx*(x[i+1][j] - 2.0*x[i][j] + x[i-1][j]);    }  }  ierr = DMDAVecRestoreArrayDOFRead(dm,Xlocal,&x);CHKERRQ(ierr);  ierr = DMDAVecRestoreArrayDOF(dm,F,&f);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm,&Xlocal);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *user){  DMDALocalInfo  info;  Field          **u,**fu;  PetscErrorCode ierr;  Vec            localX;  DM             da;  PetscFunctionBegin;  ierr = SNESGetDM(snes,(DM*)&da);CHKERRQ(ierr);  ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr);  /*     Scatter ghost points to local vector, using the 2-step process        DMGlobalToLocalBegin(), DMGlobalToLocalEnd().  */  ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,localX,&u);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,F,&fu);CHKERRQ(ierr);  ierr = FormFunctionLocal(&info,u,fu,user);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,localX,&u);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,F,&fu);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da,&localX);CHKERRQ(ierr);  PetscFunctionReturn(0); }

static PetscErrorCode SNESComputePicardJacobian_DMDA(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *mstr,void *ctx){  PetscErrorCode ierr;  DM             dm;  DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;  DMDALocalInfo  info;  Vec            Xloc;  void           *x;  PetscFunctionBegin;  if (!dmdasnes->jacobianplocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");  ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);  ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);  CHKMEMQ;  ierr = (*dmdasnes->jacobianplocal)(&info,x,*A,*B,mstr,dmdasnes->picardlocalctx);CHKERRQ(ierr);  CHKMEMQ;  ierr  = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);  ierr  = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);  *mstr = SAME_NONZERO_PATTERN;  if (*A != *B) {    ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);    ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  }  PetscFunctionReturn(0);}

/* Caculate the source term of the equations, which includes all terms except convetion term, diffusion term, and timedependent term.*/PetscErrorCode CaculateLocalSourceTerm(DM dm, Vec locX, Vec F, User user){  PetscErrorCode    ierr;  DM                dmGrad = user->dmGrad;  const PetscScalar *x;  PetscScalar       *f;  PetscInt          cStart, cell;  const PetscScalar *cellgeom;  const CellGeom    *cg;  Vec               locGrad, Grad;  const PetscScalar *grad;  DM                dmCell;  PetscFunctionBeginUser;  ierr = VecGetDM(user->cellgeom,&dmCell);CHKERRQ(ierr);  ierr = VecGetArrayRead(locX,&x);CHKERRQ(ierr);  ierr = VecGetArray(F,&f);CHKERRQ(ierr);  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr);  ierr = VecGetArrayRead(user->cellgeom,&cellgeom);CHKERRQ(ierr);  ierr = DMGetGlobalVector(dmGrad,&Grad);CHKERRQ(ierr);  ierr = DMGetLocalVector(dmGrad,&locGrad);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(dmGrad,Grad,INSERT_VALUES,locGrad);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(dmGrad,Grad,INSERT_VALUES,locGrad);CHKERRQ(ierr);  ierr = DMRestoreGlobalVector(dmGrad,&Grad);CHKERRQ(ierr);  ierr = VecGetArrayRead(locGrad,&grad);CHKERRQ(ierr);  for (cell = cStart; cell < user->cEndInterior; cell++) {    PetscScalar       *fref;    const PetscScalar *xref;    PetscScalar       *cgrad;    ierr = DMPlexPointLocalRead(dmCell,cell,cellgeom,&cg);CHKERRQ(ierr);    ierr = DMPlexPointLocalRead(dm,cell,x,&xref);CHKERRQ(ierr); /*For the unkown variables*/    ierr = DMPlexPointGlobalRef(dm,cell,f,&fref);CHKERRQ(ierr);    ierr = DMPlexPointLocalRead(dmGrad,cell,grad,&cgrad);CHKERRQ(ierr);//    if (!fref){ PetscPrintf(PETSC_COMM_WORLD,"%d, %d/n", cell, user->cEndInterior);}    if (fref){      fref[0] += SourceRho(user, cgrad, xref, cg->centroid);/*the continuity equation*/      fref[1] += SourceU(user, cgrad, xref, cg->centroid); /*Momentum U*/      fref[2] += SourceV(user, cgrad, xref, cg->centroid); /*Momentum V*/      fref[3] += SourceW(user, cgrad, xref, cg->centroid); /*Momentum W*/      fref[4] += SourceE(user, cgrad, xref, cg->centroid);/*Energy*/    }  }  ierr = VecRestoreArrayRead(locX,&x);CHKERRQ(ierr);  ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);  ierr = VecRestoreArrayRead(user->cellgeom,&cellgeom);CHKERRQ(ierr);  ierr = VecRestoreArrayRead(locGrad,&grad);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dmGrad,&locGrad);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode SetCoordinates1d(DM da){  PetscErrorCode ierr;  PetscInt       i,start,m;  Vec            local,global;  PetscScalar    *coors,*coorslocal;  DM             cda;  PetscFunctionBeginUser;  ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr);  ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr);  ierr = DMGetCoordinates(da,&global);CHKERRQ(ierr);  ierr = DMGetCoordinatesLocal(da,&local);CHKERRQ(ierr);  ierr = DMDAVecGetArray(cda,global,&coors);CHKERRQ(ierr);  ierr = DMDAVecGetArrayRead(cda,local,&coorslocal);CHKERRQ(ierr);  ierr = DMDAGetCorners(cda,&start,0,0,&m,0,0);CHKERRQ(ierr);  for (i=start; i<start+m; i++) {    if (i % 2) {      coors[i] = coorslocal[i-1] + .1*(coorslocal[i+1] - coorslocal[i-1]);    }  }  ierr = DMDAVecRestoreArray(cda,global,&coors);CHKERRQ(ierr);  ierr = DMDAVecRestoreArrayRead(cda,local,&coorslocal);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(cda,global,INSERT_VALUES,local);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(cda,global,INSERT_VALUES,local);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode SetCoordinates2d(DM da){  PetscErrorCode ierr;  PetscInt       i,j,mstart,m,nstart,n;  Vec            local,global;  DMDACoor2d     **coors,**coorslocal;  DM             cda;  PetscFunctionBeginUser;  ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr);  ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr);  ierr = DMGetCoordinates(da,&global);CHKERRQ(ierr);  ierr = DMGetCoordinatesLocal(da,&local);CHKERRQ(ierr);  ierr = DMDAVecGetArray(cda,global,&coors);CHKERRQ(ierr);  ierr = DMDAVecGetArrayRead(cda,local,&coorslocal);CHKERRQ(ierr);  ierr = DMDAGetCorners(cda,&mstart,&nstart,0,&m,&n,0);CHKERRQ(ierr);  for (i=mstart; i<mstart+m; i++) {    for (j=nstart; j<nstart+n; j++) {      if (i % 2) {        coors[j][i].x = coorslocal[j][i-1].x + .1*(coorslocal[j][i+1].x - coorslocal[j][i-1].x);      }      if (j % 2) {        coors[j][i].y = coorslocal[j-1][i].y + .3*(coorslocal[j+1][i].y - coorslocal[j-1][i].y);      }    }  }  ierr = DMDAVecRestoreArray(cda,global,&coors);CHKERRQ(ierr);  ierr = DMDAVecRestoreArrayRead(cda,local,&coorslocal);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(cda,global,INSERT_VALUES,local);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(cda,global,INSERT_VALUES,local);CHKERRQ(ierr);  PetscFunctionReturn(0);}

static PetscErrorCode TSComputeRHSJacobian_DMDA(TS ts,PetscReal ptime,Vec X,Mat *A,Mat *B,MatStructure *mstr,void *ctx){  PetscErrorCode ierr;  DM             dm;  DMTS_DA        *dmdats = (DMTS_DA*)ctx;  DMDALocalInfo  info;  Vec            Xloc;  void           *x;  PetscFunctionBegin;  if (!dmdats->rhsfunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context");  ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);  if (dmdats->rhsjacobianlocal) {    ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);    ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);    ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);    ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);    CHKMEMQ;    ierr = (*dmdats->rhsjacobianlocal)(&info,ptime,x,*A,*B,mstr,dmdats->rhsjacobianlocalctx);CHKERRQ(ierr);    CHKMEMQ;    ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);    ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);  } else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"TSComputeRHSJacobian_DMDA() called without calling DMDATSSetRHSJacobian()");  /* This will be redundant if the user called both, but it's too common to forget. */  if (*A != *B) {    ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);    ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  }  PetscFunctionReturn(0);}

/*   ComputeFunction - Evaluates nonlinear function, F(x).   Input Parameters:.  X - input vector.  user - user-defined application context   Output Parameter:.  F - function vector */PetscErrorCode ComputeFunction(AppCtx *user,Vec X,Vec F){  PetscErrorCode ierr;  PetscInt       i,j,row,mx,my,xs,ys,xm,ym,gxs,gys,gxm,gym;  PetscReal      two = 2.0,one = 1.0,lambda,hx,hy,hxdhy,hydhx,sc;  PetscScalar    u,uxx,uyy,*x,*f;  Vec            localX = user->localX;  mx = user->mx;            my = user->my;            lambda = user->param;  hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);  sc = hx*hy*lambda;        hxdhy = hx/hy;            hydhx = hy/hx;  /*     Scatter ghost points to local vector, using the 2-step process        DMGlobalToLocalBegin(), DMGlobalToLocalEnd().     By placing code between these two statements, computations can be     done while messages are in transition.  */  ierr = DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  /*     Get pointers to vector data  */  ierr = VecGetArray(localX,&x);CHKERRQ(ierr);  ierr = VecGetArray(F,&f);CHKERRQ(ierr);  /*     Get local grid boundaries  */  ierr = DMDAGetCorners(user->da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);  ierr = DMDAGetGhostCorners(user->da,&gxs,&gys,NULL,&gxm,&gym,NULL);CHKERRQ(ierr);  /*     Compute function over the locally owned part of the grid  */  for (j=ys; j<ys+ym; j++) {    row = (j - gys)*gxm + xs - gxs - 1;    for (i=xs; i<xs+xm; i++) {      row++;      if (i == 0 || j == 0 || i == mx-1 || j == my-1) {        f[row] = x[row];        continue;      }      u      = x[row];      uxx    = (two*u - x[row-1] - x[row+1])*hydhx;      uyy    = (two*u - x[row-gxm] - x[row+gxm])*hxdhy;      f[row] = uxx + uyy - sc*PetscExpScalar(u);    }  }  /*     Restore vectors  */  ierr = VecRestoreArray(localX,&x);CHKERRQ(ierr);  ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);  ierr = PetscLogFlops(11.0*ym*xm);CHKERRQ(ierr);  return 0;}

static PetscErrorCode SNESComputeJacobian_DMLocal(SNES snes,Vec X,Mat A,Mat B,void *ctx){  PetscErrorCode ierr;  DM             dm;  DMSNES_Local   *dmlocalsnes = (DMSNES_Local*)ctx;  Vec            Xloc;  PetscFunctionBegin;  ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);  if (dmlocalsnes->jacobianlocal) {    ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);    ierr = VecZeroEntries(Xloc);CHKERRQ(ierr);    if (dmlocalsnes->boundarylocal) {ierr = (*dmlocalsnes->boundarylocal)(dm,Xloc,dmlocalsnes->boundarylocalctx);CHKERRQ(ierr);}    ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);    ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);    CHKMEMQ;    ierr = (*dmlocalsnes->jacobianlocal)(dm,Xloc,A,B,dmlocalsnes->jacobianlocalctx);CHKERRQ(ierr);    CHKMEMQ;    ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);  } else {    MatFDColoring fdcoloring;    ierr = PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring);CHKERRQ(ierr);    if (!fdcoloring) {      ISColoring coloring;      ierr = DMCreateColoring(dm,dm->coloringtype,&coloring);CHKERRQ(ierr);      ierr = MatFDColoringCreate(B,coloring,&fdcoloring);CHKERRQ(ierr);      ierr = ISColoringDestroy(&coloring);CHKERRQ(ierr);      switch (dm->coloringtype) {      case IS_COLORING_GLOBAL:        ierr = MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))SNESComputeFunction_DMLocal,dmlocalsnes);CHKERRQ(ierr);        break;      default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]);      }      ierr = PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix);CHKERRQ(ierr);      ierr = MatFDColoringSetFromOptions(fdcoloring);CHKERRQ(ierr);      ierr = MatFDColoringSetUp(B,coloring,fdcoloring);CHKERRQ(ierr);      ierr = PetscObjectCompose((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject)fdcoloring);CHKERRQ(ierr);      ierr = PetscObjectDereference((PetscObject)fdcoloring);CHKERRQ(ierr);      /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call       * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the       * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually       * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be       * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.       */      ierr = PetscObjectDereference((PetscObject)dm);CHKERRQ(ierr);    }    ierr  = MatFDColoringApply(B,fdcoloring,X,snes);CHKERRQ(ierr);  }  /* This will be redundant if the user called both, but it's too common to forget. */  if (A != B) {    ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);    ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  }  PetscFunctionReturn(0);}

PetscErrorCode DMGlobalToLocalBegin_Network(DM dm, Vec g, InsertMode mode, Vec l){  PetscErrorCode ierr;  DM_Network     *network = (DM_Network*) dm->data;  PetscFunctionBegin;  ierr = DMGlobalToLocalBegin(network->plex,g,mode,l);CHKERRQ(ierr);  PetscFunctionReturn(0);}

static PetscErrorCode SNESComputeFunction_DMDA(SNES snes,Vec X,Vec F,void *ctx){  PetscErrorCode ierr;  DM             dm;  DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;  DMDALocalInfo  info;  Vec            Xloc;  void           *x,*f;  PetscFunctionBegin;  PetscValidHeaderSpecific(snes,SNES_CLASSID,1);  PetscValidHeaderSpecific(X,VEC_CLASSID,2);  PetscValidHeaderSpecific(F,VEC_CLASSID,3);  if (!dmdasnes->residuallocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");  ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);  ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);  switch (dmdasnes->residuallocalimode) {  case INSERT_VALUES: {    ierr = DMDAVecGetArray(dm,F,&f);CHKERRQ(ierr);    ierr = PetscLogEventBegin(SNES_FunctionEval,snes,X,F,0);CHKERRQ(ierr);    CHKMEMQ;    ierr = (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);CHKERRQ(ierr);    CHKMEMQ;    ierr = PetscLogEventEnd(SNES_FunctionEval,snes,X,F,0);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dm,F,&f);CHKERRQ(ierr);  } break;  case ADD_VALUES: {    Vec Floc;    ierr = DMGetLocalVector(dm,&Floc);CHKERRQ(ierr);    ierr = VecZeroEntries(Floc);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dm,Floc,&f);CHKERRQ(ierr);    ierr = PetscLogEventBegin(SNES_FunctionEval,snes,X,F,0);CHKERRQ(ierr);    CHKMEMQ;    ierr = (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);CHKERRQ(ierr);    CHKMEMQ;    ierr = PetscLogEventEnd(SNES_FunctionEval,snes,X,F,0);CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(dm,Floc,&f);CHKERRQ(ierr);    ierr = VecZeroEntries(F);CHKERRQ(ierr);    ierr = DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);    ierr = DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);    ierr = DMRestoreLocalVector(dm,&Floc);CHKERRQ(ierr);  } break;  default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode);  }  ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);  if (snes->domainerror) {    ierr = VecSetInf(F);CHKERRQ(ierr);  }  PetscFunctionReturn(0);}

/*   FormFunction - Evaluates nonlinear function, F(x).   Input Parameters:.  ts - the TS context.  X - input vector.  ptr - optional user-defined context, as set by SNESSetFunction()   Output Parameter:.  F - function vector */PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr){  DM             da;  PetscErrorCode ierr;  PetscInt       i,Mx,xs,xm;  PetscReal      hx,sx;  PetscScalar    *x,*f;  Vec            localX;  UserCtx        *ctx = (UserCtx*)ptr;  PetscFunctionBegin;  ierr = TSGetDM(ts,&da);CHKERRQ(ierr);  ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr);  ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr);  hx = 1.0/(PetscReal)Mx; sx = 1.0/(hx*hx);  /*     Scatter ghost points to local vector,using the 2-step process        DMGlobalToLocalBegin(),DMGlobalToLocalEnd().     By placing code between these two statements, computations can be     done while messages are in transition.  */  ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  /*     Get pointers to vector data  */  ierr = DMDAVecGetArrayRead(da,localX,&x);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);  /*     Get local grid boundaries  */  ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);  /*     Compute function over the locally owned part of the grid  */  for (i=xs; i<xs+xm; i++) {    f[i] = ctx->kappa*(x[i-1] + x[i+1] - 2.0*x[i])*sx;    if (ctx->allencahn) f[i] += (x[i] - x[i]*x[i]*x[i]);  }  /*     Restore vectors  */  ierr = DMDAVecRestoreArrayRead(da,localX,&x);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da,&localX);CHKERRQ(ierr);  PetscFunctionReturn(0);}

static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr){  User           user = (User)ptr;  DM             da;  DMDALocalInfo  info;  PetscInt       i;  Field          *x,*xdot,*f;  PetscReal      hx;  Vec            Xloc;  PetscErrorCode ierr;  PetscFunctionBeginUser;  ierr = TSGetDM(ts,&da);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);  hx   = 1.0/(PetscReal)(info.mx-1);  /*     Scatter ghost points to local vector,using the 2-step process        DMGlobalToLocalBegin(),DMGlobalToLocalEnd().     By placing code between these two statements, computations can be     done while messages are in transition.  */  ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  /* Get pointers to vector data */  ierr = DMDAVecGetArray(da,Xloc,&x);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,Xdot,&xdot);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);  /* Compute function over the locally owned part of the grid */  for (i=info.xs; i<info.xs+info.xm; i++) {    if (i == 0) {      f[i].u = hx * (x[i].u - user->uleft);      f[i].v = hx * (x[i].v - user->vleft);    } else if (i == info.mx-1) {      f[i].u = hx * (x[i].u - user->uright);      f[i].v = hx * (x[i].v - user->vright);    } else {      f[i].u = hx * xdot[i].u - user->alpha * (x[i-1].u - 2.*x[i].u + x[i+1].u) / hx;      f[i].v = hx * xdot[i].v - user->alpha * (x[i-1].v - 2.*x[i].v + x[i+1].v) / hx;    }  }  /* Restore vectors */  ierr = DMDAVecRestoreArray(da,Xloc,&x);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,Xdot,&xdot);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr);  PetscFunctionReturn(0);}

static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr){  User           user = (User)ptr;  DM             da;  Vec            Xloc;  DMDALocalInfo  info;  PetscInt       i,j;  PetscReal      hx;  Field          *f;  const Field    *x;  PetscErrorCode ierr;  PetscFunctionBeginUser;  ierr = TSGetDM(ts,&da);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);  hx   = 1.0/(PetscReal)info.mx;  /*     Scatter ghost points to local vector,using the 2-step process        DMGlobalToLocalBegin(),DMGlobalToLocalEnd().     By placing code between these two statements, computations can be     done while messages are in transition.  */  ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  /* Get pointers to vector data */  ierr = DMDAVecGetArrayRead(da,Xloc,(void*)&x);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);  /* Compute function over the locally owned part of the grid */  for (i=info.xs; i<info.xs+info.xm; i++) {    const PetscReal *a  = user->a;    PetscReal       u0t[2];    u0t[0] = 1.0 - PetscPowRealInt(PetscSinReal(12*t),4);    u0t[1] = 0.0;    for (j=0; j<2; j++) {      if (i == 0)              f[i][j] = a[j]/hx*(1./3*u0t[j] + 0.5*x[i][j] - x[i+1][j] + 1./6*x[i+2][j]);      else if (i == 1)         f[i][j] = a[j]/hx*(-1./12*u0t[j] + 2./3*x[i-1][j] - 2./3*x[i+1][j] + 1./12*x[i+2][j]);      else if (i == info.mx-2) f[i][j] = a[j]/hx*(-1./6*x[i-2][j] + x[i-1][j] - 0.5*x[i][j] - 1./3*x[i+1][j]);      else if (i == info.mx-1) f[i][j] = a[j]/hx*(-x[i][j] + x[i-1][j]);      else                     f[i][j] = a[j]/hx*(-1./12*x[i-2][j] + 2./3*x[i-1][j] - 2./3*x[i+1][j] + 1./12*x[i+2][j]);    }  }  /* Restore vectors */  ierr = DMDAVecRestoreArrayRead(da,Xloc,(void*)&x);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr);  PetscFunctionReturn(0);}

/*  This function should eventually replace:    DMDAComputeFunction() and DMDAComputeFunction1() */static PetscErrorCode TSComputeIFunction_DMDA(TS ts,PetscReal ptime,Vec X,Vec Xdot,Vec F,void *ctx){  PetscErrorCode ierr;  DM             dm;  DMTS_DA        *dmdats = (DMTS_DA*)ctx;  DMDALocalInfo  info;  Vec            Xloc;  void           *x,*f,*xdot;  PetscFunctionBegin;  PetscValidHeaderSpecific(ts,TS_CLASSID,1);  PetscValidHeaderSpecific(X,VEC_CLASSID,2);  PetscValidHeaderSpecific(F,VEC_CLASSID,3);  if (!dmdats->ifunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context");  ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);  ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);  ierr = DMDAVecGetArray(dm,Xdot,&xdot);CHKERRQ(ierr);  switch (dmdats->ifunctionlocalimode) {  case INSERT_VALUES: {    ierr = DMDAVecGetArray(dm,F,&f);CHKERRQ(ierr);    CHKMEMQ;    ierr = (*dmdats->ifunctionlocal)(&info,ptime,x,xdot,f,dmdats->ifunctionlocalctx);CHKERRQ(ierr);    CHKMEMQ;    ierr = DMDAVecRestoreArray(dm,F,&f);CHKERRQ(ierr);  } break;  case ADD_VALUES: {    Vec Floc;    ierr = DMGetLocalVector(dm,&Floc);CHKERRQ(ierr);    ierr = VecZeroEntries(Floc);CHKERRQ(ierr);    ierr = DMDAVecGetArray(dm,Floc,&f);CHKERRQ(ierr);    CHKMEMQ;    ierr = (*dmdats->ifunctionlocal)(&info,ptime,x,xdot,f,dmdats->ifunctionlocalctx);CHKERRQ(ierr);    CHKMEMQ;    ierr = DMDAVecRestoreArray(dm,Floc,&f);CHKERRQ(ierr);    ierr = VecZeroEntries(F);CHKERRQ(ierr);    ierr = DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);    ierr = DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);    ierr = DMRestoreLocalVector(dm,&Floc);CHKERRQ(ierr);  } break;  default: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdats->ifunctionlocalimode);  }  ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(dm,Xdot,&xdot);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode SetInitialValues(DM networkdm,Vec X,void* appctx){  PetscErrorCode ierr;  VERTEXDATA     bus;  GEN            gen;  PetscInt       v, vStart, vEnd, offset;  PetscBool      ghostvtex;  Vec            localX;  PetscScalar    *xarr;  PetscInt       key,numComps,j,offsetd;  DMNetworkComponentGenericDataType *arr;    PetscFunctionBegin;  ierr = DMNetworkGetVertexRange(networkdm,&vStart, &vEnd);CHKERRQ(ierr);  ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr);  ierr = VecSet(X,0.0);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = VecGetArray(localX,&xarr);CHKERRQ(ierr);  ierr = DMNetworkGetComponentDataArray(networkdm,&arr);CHKERRQ(ierr);  for (v = vStart; v < vEnd; v++) {    ierr = DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);CHKERRQ(ierr);    if (ghostvtex) continue;       ierr = DMNetworkGetVariableOffset(networkdm,v,&offset);CHKERRQ(ierr);    ierr = DMNetworkGetNumComponents(networkdm,v,&numComps);CHKERRQ(ierr);    for (j=0; j < numComps; j++) {      ierr = DMNetworkGetComponentTypeOffset(networkdm,v,j,&key,&offsetd);CHKERRQ(ierr);      if (key == 1) {	bus = (VERTEXDATA)(arr+offsetd);	xarr[offset] = bus->va*PETSC_PI/180.0;	xarr[offset+1] = bus->vm;      } else if(key == 2) {	gen = (GEN)(arr+offsetd);	if (!gen->status) continue;	xarr[offset+1] = gen->vs; 	break;      }    }  }  ierr = VecRestoreArray(localX,&xarr);CHKERRQ(ierr);  ierr = DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(networkdm,&localX);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode FormJacobian_Grid(AppCtx *user,GridCtx *grid,Vec X, Mat *J,Mat *B){  Mat            jac = *J;  PetscErrorCode ierr;  PetscInt       i, j, row, mx, my, xs, ys, xm, ym, Xs, Ys, Xm, Ym, col[5], nloc, *ltog, grow;  PetscScalar    two = 2.0, one = 1.0, lambda, v[5], hx, hy, hxdhy, hydhx, sc, *x, value;  Vec            localX = grid->localX;  mx = grid->mx;            my = grid->my;            lambda = user->param;  hx = one/(double)(mx-1);  hy = one/(double)(my-1);  sc = hx*hy;               hxdhy = hx/hy;            hydhx = hy/hx;  /* Get ghost points */  ierr = DMGlobalToLocalBegin(grid->da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(grid->da,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMDAGetCorners(grid->da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);  ierr = DMDAGetGhostCorners(grid->da,&Xs,&Ys,0,&Xm,&Ym,0);CHKERRQ(ierr);  ierr = DMDAGetGlobalIndices(grid->da,&nloc,&ltog);CHKERRQ(ierr);  ierr = VecGetArray(localX,&x);CHKERRQ(ierr);  /* Evaluate Jacobian of function */  for (j=ys; j<ys+ym; j++) {    row = (j - Ys)*Xm + xs - Xs - 1;     for (i=xs; i<xs+xm; i++) {      row++;      grow = ltog[row];      if (i > 0 && i < mx-1 && j > 0 && j < my-1) {        v[0] = -hxdhy; col[0] = ltog[row - Xm];        v[1] = -hydhx; col[1] = ltog[row - 1];        v[2] = two*(hydhx + hxdhy) - sc*lambda*exp(x[row]); col[2] = grow;        v[3] = -hydhx; col[3] = ltog[row + 1];        v[4] = -hxdhy; col[4] = ltog[row + Xm];        ierr = MatSetValues(jac,1,&grow,5,col,v,INSERT_VALUES);CHKERRQ(ierr);      } else if ((i > 0 && i < mx-1) || (j > 0 && j < my-1)){        value = .5*two*(hydhx + hxdhy);        ierr = MatSetValues(jac,1,&grow,1,&grow,&value,INSERT_VALUES);CHKERRQ(ierr);      } else {        value = .25*two*(hydhx + hxdhy);        ierr = MatSetValues(jac,1,&grow,1,&grow,&value,INSERT_VALUES);CHKERRQ(ierr);      }    }  }  ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  ierr = VecRestoreArray(localX,&x);CHKERRQ(ierr);  ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  return 0;}