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

/*     Matrix free operation of 1d Laplacian and Grad for GLL spectral elements*/PetscErrorCode MatMult_Laplacian(Mat A,Vec x,Vec y){  AppCtx            *appctx;  PetscErrorCode    ierr;  PetscReal         **temp,vv;  PetscInt          i,j,xs,xn;  Vec               xlocal,ylocal;  const PetscScalar *xl;  PetscScalar       *yl;  PetscBLASInt      _One = 1,n;  PetscScalar       _DOne = 1;    ierr = MatShellGetContext(A,&appctx);CHKERRQ(ierr);  ierr = DMGetLocalVector(appctx->da,&xlocal);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(appctx->da,x,INSERT_VALUES,xlocal);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(appctx->da,x,INSERT_VALUES,xlocal);CHKERRQ(ierr);  ierr = DMGetLocalVector(appctx->da,&ylocal);CHKERRQ(ierr);  ierr = VecSet(ylocal,0.0);CHKERRQ(ierr);  ierr = PetscGLLElementLaplacianCreate(&appctx->SEMop.gll,&temp);CHKERRQ(ierr);  for (i=0; i<appctx->param.N; i++) {    vv =-appctx->param.mu*2.0/appctx->param.Le;    for (j=0; j<appctx->param.N; j++) temp[i][j]=temp[i][j]*vv;  }  ierr = DMDAVecGetArrayRead(appctx->da,xlocal,(void*)&xl);CHKERRQ(ierr);  ierr = DMDAVecGetArray(appctx->da,ylocal,&yl);CHKERRQ(ierr);  ierr = DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);CHKERRQ(ierr);  ierr = PetscBLASIntCast(appctx->param.N,&n);CHKERRQ(ierr);  for (j=xs; j<xs+xn; j += appctx->param.N-1) {    PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&_DOne,&temp[0][0],&n,&xl[j],&_One,&_DOne,&yl[j],&_One));  }  ierr = DMDAVecRestoreArrayRead(appctx->da,xlocal,(void*)&xl);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(appctx->da,ylocal,&yl);CHKERRQ(ierr);  ierr = PetscGLLElementLaplacianDestroy(&appctx->SEMop.gll,&temp);CHKERRQ(ierr);  ierr = VecSet(y,0.0);CHKERRQ(ierr);  ierr = DMLocalToGlobalBegin(appctx->da,ylocal,ADD_VALUES,y);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(appctx->da,ylocal,ADD_VALUES,y);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(appctx->da,&xlocal);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(appctx->da,&ylocal);CHKERRQ(ierr);  ierr = VecPointwiseDivide(y,y,appctx->SEMop.mass);CHKERRQ(ierr);  return 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);}

/* no discretization is given so DMGetNumFields and PetscDSGetTotalDimension yield 0 */static PetscErrorCode test2(DM dm, AppCtx *options){  IS             cells;  Vec            locX, locX_t, locA;  PetscScalar    *u, *u_t, *a;  PetscMPIInt    rank;  PetscErrorCode ierr;  PetscFunctionBegin;  ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);  ierr = ISCreateStride(PETSC_COMM_SELF, rank ? 0 : 1, 0, 1, &cells);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm, &locX);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm, &locX_t);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm, &locA);CHKERRQ(ierr);  ierr = DMPlexGetCellFields(    dm, cells, locX, locX_t, locA, &u, &u_t, &a);CHKERRQ(ierr);  ierr = DMPlexRestoreCellFields(dm, cells, locX, locX_t, locA, &u, &u_t, &a);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm, &locX);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm, &locX_t);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm, &locA);CHKERRQ(ierr);  ierr = ISDestroy(&cells);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PETSC_INTERN PetscErrorCode DMStagSetUniformCoordinatesExplicit_1d(DM dm,PetscReal xmin,PetscReal xmax){  PetscErrorCode ierr;  DM_Stag        *stagCoord;  DM             dmCoord;  Vec            coordLocal,coord;  PetscReal      h,min;  PetscScalar    **arr;  PetscInt       ind,start,n,nExtra,s;  PetscInt       ileft,ielement;  PetscFunctionBegin;  ierr = DMGetCoordinateDM(dm, &dmCoord);CHKERRQ(ierr);  stagCoord = (DM_Stag*) dmCoord->data;  for (s=0; s<2; ++s) {    if (stagCoord->dof[s] !=0 && stagCoord->dof[s] != 1) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Coordinate DM in 1 dimensions must have 0 or 1 dof on each stratum, but stratum %d has %d dof",s,stagCoord->dof[s]);  }  ierr = DMGetLocalVector(dmCoord,&coordLocal);CHKERRQ(ierr);  ierr = DMStagVecGetArrayDOF(dmCoord,coordLocal,&arr);CHKERRQ(ierr);  if (stagCoord->dof[0]) {    ierr = DMStagGetLocationSlot(dmCoord,DMSTAG_LEFT,0,&ileft);CHKERRQ(ierr);  }  if (stagCoord->dof[1]) {    ierr = DMStagGetLocationSlot(dmCoord,DMSTAG_ELEMENT,0,&ielement);CHKERRQ(ierr);  }  ierr = DMStagGetCorners(dmCoord,&start,NULL,NULL,&n,NULL,NULL,&nExtra,NULL,NULL);CHKERRQ(ierr);  min = xmin;  h = (xmax-xmin)/stagCoord->N[0];  for(ind=start; ind<start + n + nExtra; ++ind) {    if (stagCoord->dof[0]) {      const PetscReal off = 0.0;        arr[ind][ileft] = min + ((PetscReal)ind + off) * h;    }    if (stagCoord->dof[1]) {      const PetscReal off = 0.5;        arr[ind][ielement] = min + ((PetscReal)ind + off) * h;    }  }  ierr = DMStagVecRestoreArrayDOF(dmCoord,coordLocal,&arr);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(dmCoord,&coord);CHKERRQ(ierr);  ierr = DMLocalToGlobalBegin(dmCoord,coordLocal,INSERT_VALUES,coord);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(dmCoord,coordLocal,INSERT_VALUES,coord);CHKERRQ(ierr);  ierr = DMSetCoordinates(dm,coord);CHKERRQ(ierr);  ierr = PetscLogObjectParent((PetscObject)dm,(PetscObject)coord);CHKERRQ(ierr);  ierr = VecDestroy(&coord);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dmCoord,&coordLocal);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode ComputePredictor(DM da, UserContext *user){  Vec            uOldLocal, uLocal,uOld;  PetscScalar    *pOld;  PetscScalar    *p;  PetscErrorCode ierr;  PetscFunctionBeginUser;  ierr = DMGetGlobalVector(da, &uOld);CHKERRQ(ierr);  ierr = DMGetLocalVector(da, &uOldLocal);CHKERRQ(ierr);  ierr = DMGetLocalVector(da, &uLocal);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(da, uOld, INSERT_VALUES, uOldLocal);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da, uOld, INSERT_VALUES, uOldLocal);CHKERRQ(ierr);  ierr = VecGetArray(uOldLocal, &pOld);CHKERRQ(ierr);  ierr = VecGetArray(uLocal,    &p);CHKERRQ(ierr);  /* Source terms are all zero right now */  ierr = CalculateElementVelocity(da, user);CHKERRQ(ierr);  ierr = TaylorGalerkinStepI(da, user);CHKERRQ(ierr);  ierr = TaylorGalerkinStepIIMomentum(da, user);CHKERRQ(ierr);  ierr = TaylorGalerkinStepIIMassEnergy(da, user);CHKERRQ(ierr);  /* Solve equation (9) for /delta(/rho/vu) and (/rho/vu)^* */  /* Solve equation (13) for /delta/rho and /rho^* */  /* Solve equation (15) for /delta(/rho e_t) and (/rho e_t)^* */  /* Apply artifical dissipation */  /* Determine the smoothed explicit pressure, /tilde P and temperature /tilde T using the equation of state */  ierr = VecRestoreArray(uOldLocal, &pOld);CHKERRQ(ierr);  ierr = VecRestoreArray(uLocal,    &p);CHKERRQ(ierr);#if 0  ierr = DMLocalToGlobalBegin(da, uLocal, ADD_VALUES,u);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(da, uLocal, ADD_VALUES,u);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da, &uOldLocal);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da, &uLocal);CHKERRQ(ierr);#endif  PetscFunctionReturn(0);}

PetscErrorCode CreatePressureNullSpace(DM dm, AppCtx *user, MatNullSpace *nullSpace){  Vec            vec, localVec;  PetscErrorCode ierr;  PetscFunctionBeginUser;  ierr = DMGetGlobalVector(dm, &vec);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm, &localVec);CHKERRQ(ierr);  ierr = VecSet(vec,  0.0);CHKERRQ(ierr);  /* Put a constant in for all pressures     Could change this to project the constant function onto the pressure space (when that is finished) */  {    PetscSection section;    PetscInt     pStart, pEnd, p;    PetscScalar  *a;    ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);    ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);    ierr = VecGetArray(localVec, &a);CHKERRQ(ierr);    for (p = pStart; p < pEnd; ++p) {      PetscInt fDim, off, d;      ierr = PetscSectionGetFieldDof(section, p, 1, &fDim);CHKERRQ(ierr);      ierr = PetscSectionGetFieldOffset(section, p, 1, &off);CHKERRQ(ierr);      for (d = 0; d < fDim; ++d) a[off+d] = 1.0;    }    ierr = VecRestoreArray(localVec, &a);CHKERRQ(ierr);  }  ierr = DMLocalToGlobalBegin(dm, localVec, INSERT_VALUES, vec);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(dm, localVec, INSERT_VALUES, vec);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm, &localVec);CHKERRQ(ierr);  ierr = VecNormalize(vec, NULL);CHKERRQ(ierr);  if (user->debug) {    ierr = PetscPrintf(PetscObjectComm((PetscObject)dm), "Pressure Null Space/n");CHKERRQ(ierr);    ierr = VecView(vec, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);  }  ierr = MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullSpace);CHKERRQ(ierr);  ierr = DMRestoreGlobalVector(dm, &vec);CHKERRQ(ierr);  /* New style for field null spaces */  {    PetscObject  pressure;    MatNullSpace nullSpacePres;    ierr = DMGetField(dm, 1, &pressure);CHKERRQ(ierr);    ierr = MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nullSpacePres);CHKERRQ(ierr);    ierr = PetscObjectCompose(pressure, "nullspace", (PetscObject) nullSpacePres);CHKERRQ(ierr);    ierr = MatNullSpaceDestroy(&nullSpacePres);CHKERRQ(ierr);  }  PetscFunctionReturn(0);}

PetscErrorCode TestVecClosure(DM dm, AppCtx *user){  PetscSection       s;  Vec                v;  PetscInt           numRuns, cStart, cEnd, c, i;  PetscScalar        tmpArray[64];  PetscScalar       *userArray     = user->reuseArray ? tmpArray : NULL;  PetscReal          maxTimePerRun = user->maxVecClosureTime;  PetscStageLog      stageLog;  PetscEventPerfLog  eventLog;  PetscInt           stage;  PetscLogEvent      event;  PetscEventPerfInfo eventInfo;  PetscErrorCode     ierr;  PetscFunctionBegin;  ierr = PetscLogStageRegister("DMPlex Vector Closure Test", &stage);CHKERRQ(ierr);  ierr = PetscLogEventRegister("VecClosure", PETSC_OBJECT_CLASSID, &event);CHKERRQ(ierr);  ierr = PetscLogStagePush(stage);CHKERRQ(ierr);  ierr = DMPlexCreateSection(dm, user->dim, user->numFields, user->numComponents, user->numDof, 0, NULL, NULL, &s);CHKERRQ(ierr);  ierr = DMSetDefaultSection(dm, s);CHKERRQ(ierr);  ierr = PetscSectionDestroy(&s);CHKERRQ(ierr);  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm, &v);CHKERRQ(ierr);  ierr = PetscLogEventBegin(event,0,0,0,0);CHKERRQ(ierr);  for (i = 0; i < user->iterations; ++i) {    for (c = cStart; c < cEnd; ++c) {      PetscScalar *closure     = userArray;      PetscInt     closureSize = 64;;      ierr = DMPlexVecGetClosure(dm, s, v, c, &closureSize, &closure);CHKERRQ(ierr);      if (!user->reuseArray) {ierr = DMPlexVecRestoreClosure(dm, s, v, c, &closureSize, &closure);CHKERRQ(ierr);}    }  }  ierr = PetscLogEventEnd(event,0,0,0,0);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm, &v);CHKERRQ(ierr);  ierr = PetscLogStagePop();CHKERRQ(ierr);  ierr = PetscLogGetStageLog(&stageLog);  ierr = PetscStageLogGetEventPerfLog(stageLog, stage, &eventLog);  numRuns   = (cEnd-cStart) * user->iterations;  eventInfo = eventLog->eventInfo[event];  if (eventInfo.count != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of event calls %d should be %d", eventInfo.count, 1);  if ((PetscInt) eventInfo.flops != 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of event flops %d should be %d", (PetscInt) eventInfo.flops, 0);  if (eventInfo.time > maxTimePerRun * numRuns) {    ierr = PetscPrintf(PETSC_COMM_SELF, "VecClosures: %d Average time per cone: %gs standard: %gs/n", numRuns, eventInfo.time/numRuns, maxTimePerRun);    if (user->errors) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Average time for vector closure %g > standard %g", eventInfo.time/numRuns, maxTimePerRun);  }  PetscFunctionReturn(0);}

PetscErrorCode ComputeCorrector(DM da, Vec uOld, Vec u){  Vec            uOldLocal, uLocal;  PetscScalar    *cOld;  PetscScalar    *c;  PetscInt       i,ne,nc;  const PetscInt *e;  PetscErrorCode ierr;  PetscFunctionBeginUser;  ierr = VecSet(u,0.0);CHKERRQ(ierr);  ierr = DMGetLocalVector(da, &uOldLocal);CHKERRQ(ierr);  ierr = DMGetLocalVector(da, &uLocal);CHKERRQ(ierr);  ierr = VecSet(uLocal,0.0);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(da, uOld, INSERT_VALUES, uOldLocal);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da, uOld, INSERT_VALUES, uOldLocal);CHKERRQ(ierr);  ierr = VecGetArray(uOldLocal, &cOld);CHKERRQ(ierr);  ierr = VecGetArray(uLocal,    &c);CHKERRQ(ierr);  /* access the list of elements on this processor and loop over them */  ierr = DMDAGetElements(da,&ne,&nc,&e);CHKERRQ(ierr);  for (i=0; i<ne; i++) {    /* this is nonsense, but copy each nodal value*/    c[e[3*i]]   = cOld[e[3*i]];    c[e[3*i+1]] = cOld[e[3*i+1]];    c[e[3*i+2]] = cOld[e[3*i+2]];  }  ierr = DMDARestoreElements(da,&ne,&nc,&e);CHKERRQ(ierr);  ierr = VecRestoreArray(uOldLocal, &cOld);CHKERRQ(ierr);  ierr = VecRestoreArray(uLocal,    &c);CHKERRQ(ierr);  ierr = DMLocalToGlobalBegin(da, uLocal, ADD_VALUES,u);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(da, uLocal, ADD_VALUES,u);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da, &uOldLocal);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da, &uLocal);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 WaterFormFunction(SNES snes,Vec X, Vec F, void *user){  PetscErrorCode ierr;  DM             networkdm;  Vec            localX,localF;  const PetscInt *v,*e;  PetscInt       nv,ne;  PetscFunctionBegin;  /* Get the DM attached with the SNES */  ierr = SNESGetDM(snes,&networkdm);CHKERRQ(ierr);  /* Get local vertices and edges */  ierr = DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&v,&e);CHKERRQ(ierr);  /* Get local vectors */  ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr);  ierr = DMGetLocalVector(networkdm,&localF);CHKERRQ(ierr);  /* Scatter values from global to local vector */  ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);  /* Initialize residual */  ierr = VecSet(F,0.0);CHKERRQ(ierr);  ierr = VecSet(localF,0.0);CHKERRQ(ierr);  ierr = FormFunction_Water(networkdm,localX,localF,nv,ne,v,e,NULL);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(networkdm,&localX);CHKERRQ(ierr);  ierr = DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(networkdm,&localF);CHKERRQ(ierr);  PetscFunctionReturn(0);}

/* Use Lax-Wendroff method to evaluate F(u,t) = du/dt + a *  du/dx*/PetscErrorCode IFunction_LaxWendroff(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void* ctx){  PetscErrorCode ierr;  AppCtx         *appctx=(AppCtx*)ctx;  PetscInt       mstart,mend,M,i,um;  DM             da;  Vec            Uold,localUold;  PetscScalar    *uarray,*f,*uoldarray,h,RFlux,LFlux,lambda;  PetscReal      dt,a;  PetscFunctionBegin;  ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr);  ierr = TSGetSolution(ts,&Uold);CHKERRQ(ierr);  ierr = TSGetDM(ts,&da);CHKERRQ(ierr);  ierr = DMDAGetInfo(da,0,&M,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);  ierr = DMDAGetCorners(da,&mstart,0,0,&um,0,0);CHKERRQ(ierr);  h    = 1.0/M;  mend = mstart + um;  /* printf(" mstart %d, um %d/n",mstart,um); */  ierr = DMGetLocalVector(da,&localUold);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(da,Uold,INSERT_VALUES,localUold);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,Uold,INSERT_VALUES,localUold);CHKERRQ(ierr);  /* Get pointers to vector data */  ierr = DMDAVecGetArrayRead(da,U,&uarray);CHKERRQ(ierr);  ierr = DMDAVecGetArrayRead(da,localUold,&uoldarray);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);  /* advection -- finite volume (appctx->a < 0 -- can be relaxed?) */  lambda = dt/h;  a = appctx->a;  for (i=mstart; i<mend; i++) {    RFlux = 0.5 * a * (uoldarray[i+1] + uoldarray[i]) - a*a*0.5*lambda * (uoldarray[i+1] - uoldarray[i]);    LFlux = 0.5 * a * (uoldarray[i-1] + uoldarray[i]) - a*a*0.5*lambda * (uoldarray[i] - uoldarray[i-1]);    f[i]  = uarray[i] - uoldarray[i] + lambda * (RFlux - LFlux);  }  /* Restore vectors */  ierr = DMDAVecRestoreArrayRead(da,U,&uarray);CHKERRQ(ierr);  ierr = DMDAVecRestoreArrayRead(da,localUold,&uoldarray);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da,&localUold);CHKERRQ(ierr);  PetscFunctionReturn(0);}

/* Compute pointwise residual f(x) over the locally-owned part of the grid   This is a finite difference method. In TeX, the formula we compute isf_i = /eta_{i+1/2} H_{i+1/2} (u_{i+1}-u_i) - /eta_{i-1/2} H_{i-1/2} (u_i-u_{i-1})      - dx K (H_{i+1/2}^2 - H_{i-1/2}^2)where/eta_{i+1/2} = /left|/frac{u_{i+1}-u_i}{dx}/right|^{p-2}with some regularization using user.epsilon, and  dx = L / Mx  p = 1+1/n  K = rho * g * (1-r) / (4 * B)  r = rho / rhow  B = A^{1/n}*/PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar *u,PetscScalar *f,AppCtx *user){  PetscErrorCode ierr;  PetscReal      hx, p, K, B, duL, duR, sL, sR;  PetscScalar    *H;  PetscInt       i, Mx;  Vec            localH;  PetscFunctionBegin;  ierr = DMGetLocalVector(info->da,&localH);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(info->da,user->H,INSERT_VALUES,localH); CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(info->da,user->H,INSERT_VALUES,localH); CHKERRQ(ierr);  p  = 1.0 + 1.0 / user->n;  B  = PetscPowScalar(user->A,-1.0/user->n);  K  = user->rho * user->g * (1.0 - user->rho/user->rhow) / (4.0 * B);  Mx = info->mx;  hx = user->L / ((PetscReal)Mx - 1.0);  ierr = DMDAVecGetArray(info->da,localH,&H);CHKERRQ(ierr);  for (i=info->xs; i<info->xs+info->xm; i++) {    if (i == 0) {      f[0] = u[0] - user->ug;  /* Dirichlet condition */    } else {      if (i == 1) { /* use Dirichlet condition as value for neighbor, which symmetrizes */        duL = u[i] - user->ug;      } else {        duL = u[i] - u[i-1];      }      sL = GetEta(duL, hx, p, user->epsilon) * H[i-1] * duL;      if (i == Mx-1) {  /* Neumann: calving front stress boundary condition */        duR = u[Mx-2] + 2.0 * hx * user->gamma - u[Mx-1];      } else {        duR = u[i+1] - u[i];      }      sR = GetEta(duR, hx, p, user->epsilon) * H[i] * duR;      f[i] = sR - sL - hx * K * (H[i]*H[i] - H[i-1]*H[i-1]);    }  }  ierr = DMDAVecRestoreArray(info->da,localH,&H);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(info->da,&localH);CHKERRQ(ierr);  PetscFunctionReturn(0);}

static PetscErrorCode test5(DM dm, AppCtx *options){  IS             cells;  Vec            locX, locX_t, locA;  PetscScalar    *u, *u_t, *a;  PetscErrorCode ierr;  PetscFunctionBegin;  locX_t = NULL;  locA = NULL;  ierr = ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &cells);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm, &locX);CHKERRQ(ierr);  ierr = DMPlexGetCellFields(    dm, cells, locX, locX_t, locA, &u, &u_t, &a);CHKERRQ(ierr);  ierr = DMPlexRestoreCellFields(dm, cells, locX, locX_t, locA, &u, &u_t, &a);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm, &locX);CHKERRQ(ierr);  ierr = ISDestroy(&cells);CHKERRQ(ierr);  PetscFunctionReturn(0);}

/*@  DMPlexTSComputeRHSFunctionFVM - Form the local forcing F from the local input X using pointwise functions specified by the user  Input Parameters:+ dm - The mesh. t - The time. locX  - Local solution- user - The user context  Output Parameter:. F  - Global output vector  Level: developer.seealso: DMPlexComputeJacobianActionFEM()@*/PetscErrorCode DMPlexTSComputeRHSFunctionFVM(DM dm, PetscReal time, Vec locX, Vec F, void *user){  Vec            locF;  PetscInt       cStart, cEnd, cEndInterior;  PetscErrorCode ierr;  PetscFunctionBegin;  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);  ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr);  cEnd = cEndInterior < 0 ? cEnd : cEndInterior;  ierr = DMGetLocalVector(dm, &locF);CHKERRQ(ierr);  ierr = VecZeroEntries(locF);CHKERRQ(ierr);  ierr = DMPlexComputeResidual_Internal(dm, cStart, cEnd, time, locX, NULL, locF, user);CHKERRQ(ierr);  ierr = DMLocalToGlobalBegin(dm, locF, INSERT_VALUES, F);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(dm, locF, INSERT_VALUES, F);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm, &locF);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode MatMult_DAAD(Mat A,Vec xx,Vec yy){    Mat_DAAD *a = (Mat_DAAD*)A->data;    Vec      localxx;    PetscErrorCode ierr;    PetscFunctionBegin;    ierr = DMGetLocalVector(a->da,&localxx);    CHKERRQ(ierr);    ierr = DMGlobalToLocalBegin(a->da,xx,INSERT_VALUES,localxx);    CHKERRQ(ierr);    ierr = DMGlobalToLocalEnd(a->da,xx,INSERT_VALUES,localxx);    CHKERRQ(ierr);    ierr = DMDAMultiplyByJacobian1WithAD(a->da,a->localu,localxx,yy,a->ctx);    CHKERRQ(ierr);    ierr = DMRestoreLocalVector(a->da,&localxx);    CHKERRQ(ierr);    PetscFunctionReturn(0);}

/**This function is for the explicit method, that is the right hand sideof the system: du/dt = F(u)*/PetscErrorCode MyRHSFunction(TS ts,PetscReal time,Vec in,Vec out,void *ctx){  User           user = (User)ctx;  DM             dm, dmFace, dmCell;  PetscSection   section;  Vec            inLocal;  PetscErrorCode ierr;  PetscFunctionBeginUser;  {    PetscReal norm;    PetscInt size;    ierr = VecNorm(in,NORM_INFINITY,&norm);CHKERRQ(ierr);    ierr = VecGetSize(in, &size);CHKERRQ(ierr);    norm = norm/size;    if (norm>1.e5) {      SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_LIB,      "The norm of the solution is: %f (current time: %f). The explicit method is going to DIVERGE!!!", norm, time);    }  }  ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);  ierr = VecGetDM(user->facegeom,&dmFace);CHKERRQ(ierr);  ierr = VecGetDM(user->cellgeom,&dmCell);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm,&inLocal);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(dm, in, INSERT_VALUES, inLocal);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(dm, in, INSERT_VALUES, inLocal);CHKERRQ(ierr);  ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);  ierr = ApplyBC(dm, time, inLocal, user);CHKERRQ(ierr);  ierr = ConstructCellCentriodGradient(user->dm, dmFace, dmCell, time, inLocal, out, user);CHKERRQ(ierr);  ierr = VecZeroEntries(out);CHKERRQ(ierr);  if (user->second_order){    ierr = CaculateLocalFunction_LS(user->dm, dmFace, dmCell, time, inLocal, out, user);CHKERRQ(ierr);  }else{    ierr = CaculateLocalFunction_Upwind(user->dm, dmFace, dmCell, time, inLocal, out, user);CHKERRQ(ierr);  }  ierr = DMRestoreLocalVector(dm,&inLocal);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode MatMult_GlobalToLocalNormal(Mat CtC, Vec x, Vec y){  DM                    dm;  Vec                   local, mask;  projectConstraintsCtx *ctx;  PetscErrorCode        ierr;  PetscFunctionBegin;  ierr = MatShellGetContext(CtC,&ctx);CHKERRQ(ierr);  dm   = ctx->dm;  mask = ctx->mask;  ierr = DMGetLocalVector(dm,&local);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(dm,x,INSERT_VALUES,local);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(dm,x,INSERT_VALUES,local);CHKERRQ(ierr);  if (mask) {ierr = VecPointwiseMult(local,mask,local);CHKERRQ(ierr);}  ierr = VecSet(y,0.);CHKERRQ(ierr);  ierr = DMLocalToGlobalBegin(dm,local,ADD_VALUES,y);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(dm,local,ADD_VALUES,y);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm,&local);CHKERRQ(ierr);  PetscFunctionReturn(0);}

/*@  DMPlexTSComputeRHSFunctionFVM - Form the local forcing F from the local input X using pointwise functions specified by the user  Input Parameters:+ dm - The mesh. t - The time. locX  - Local solution- user - The user context  Output Parameter:. F  - Global output vector  Level: developer.seealso: DMPlexComputeJacobianActionFEM()@*/PetscErrorCode DMPlexTSComputeRHSFunctionFVM(DM dm, PetscReal time, Vec locX, Vec F, void *user){  Vec            locF;  PetscInt       cStart, cEnd, cEndInterior;  DM             plex;  PetscErrorCode ierr;  PetscFunctionBegin;  ierr = DMTSConvertPlex(dm,&plex,PETSC_TRUE);CHKERRQ(ierr);  ierr = DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);CHKERRQ(ierr);  ierr = DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr);  cEnd = cEndInterior < 0 ? cEnd : cEndInterior;  ierr = DMGetLocalVector(plex, &locF);CHKERRQ(ierr);  ierr = VecZeroEntries(locF);CHKERRQ(ierr);  ierr = DMPlexComputeResidual_Internal(plex, cStart, cEnd, time, locX, NULL, time, locF, user);CHKERRQ(ierr);  ierr = DMLocalToGlobalBegin(plex, locF, ADD_VALUES, F);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(plex, locF, ADD_VALUES, F);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(plex, &locF);CHKERRQ(ierr);  ierr = DMDestroy(&plex);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode DMVecViewLocal(DM dm, Vec v, PetscViewer viewer){  Vec            lv;  PetscInt       p;  PetscMPIInt    rank, numProcs;  PetscErrorCode ierr;  PetscFunctionBeginUser;  ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);  ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm, &lv);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(dm, v, INSERT_VALUES, lv);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(dm, v, INSERT_VALUES, lv);CHKERRQ(ierr);  ierr = PetscPrintf(PETSC_COMM_WORLD, "Local function/n");CHKERRQ(ierr);  for (p = 0; p < numProcs; ++p) {    if (p == rank) {ierr = VecView(lv, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);}    ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr);  }  ierr = DMRestoreLocalVector(dm, &lv);CHKERRQ(ierr);  PetscFunctionReturn(0);}

/*@C    DMCompositeRestoreLocalVectors - Restores local vectors for each part of a DMComposite.    Not Collective    Input Parameter:.    dm - the packer object    Output Parameter:.   Vec ... - the individual sequential Vecs    Level: advanced.seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),         DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),         DMCompositeGetLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()@*/PetscErrorCode  DMCompositeRestoreLocalVectors(DM dm,...){  va_list                Argp;  PetscErrorCode         ierr;  struct DMCompositeLink *next;  DM_Composite           *com = (DM_Composite*)dm->data;  PetscFunctionBegin;  PetscValidHeaderSpecific(dm,DM_CLASSID,1);  next = com->next;  /* loop over packed objects, handling one at at time */  va_start(Argp,dm);  while (next) {    Vec *vec;    vec = va_arg(Argp, Vec*);    if (vec) {ierr = DMRestoreLocalVector(next->dm,vec);CHKERRQ(ierr);}    next = next->next;  }  va_end(Argp);  PetscFunctionReturn(0);}

static PetscErrorCode TestFunctionProjection(DM dm, DM auxdm, DMLabel label, Vec la, const char name[], AppCtx *user){  PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *);  Vec               x, lx;  PetscInt          Nf, f;  PetscInt          val[1] = {1};  char              lname[PETSC_MAX_PATH_LEN];  PetscErrorCode    ierr;  PetscFunctionBeginUser;  if (auxdm) {    ierr = PetscObjectCompose((PetscObject) dm, "dmAux", (PetscObject) auxdm);CHKERRQ(ierr);    ierr = PetscObjectCompose((PetscObject) dm, "A", (PetscObject) la);CHKERRQ(ierr);  }  ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr);  ierr = PetscMalloc1(Nf, &funcs);CHKERRQ(ierr);  for (f = 0; f < Nf; ++f) funcs[f] = linear;  ierr = DMGetGlobalVector(dm, &x);CHKERRQ(ierr);  ierr = PetscStrcpy(lname, "Function ");CHKERRQ(ierr);  ierr = PetscStrcat(lname, name);CHKERRQ(ierr);  ierr = PetscObjectSetName((PetscObject) x, lname);CHKERRQ(ierr);  if (!label) {ierr = DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_VALUES, x);CHKERRQ(ierr);}  else        {ierr = DMProjectFunctionLabel(dm, 0.0, label, 1, val, 0, NULL, funcs, NULL, INSERT_VALUES, x);CHKERRQ(ierr);}  ierr = VecViewFromOptions(x, NULL, "-func_view");CHKERRQ(ierr);  ierr = DMRestoreGlobalVector(dm, &x);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm, &lx);CHKERRQ(ierr);  ierr = PetscStrcpy(lname, "Local Function ");CHKERRQ(ierr);  ierr = PetscStrcat(lname, name);CHKERRQ(ierr);  ierr = PetscObjectSetName((PetscObject) lx, lname);CHKERRQ(ierr);  if (!label) {ierr = DMProjectFunctionLocal(dm, 0.0, funcs, NULL, INSERT_VALUES, lx);CHKERRQ(ierr);}  else        {ierr = DMProjectFunctionLabelLocal(dm, 0.0, label, 1, val, 0, NULL, funcs, NULL, INSERT_VALUES, lx);CHKERRQ(ierr);}  ierr = VecViewFromOptions(lx, NULL, "-local_func_view");CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm, &lx);CHKERRQ(ierr);  ierr = PetscFree(funcs);CHKERRQ(ierr);  if (auxdm) {    ierr = PetscObjectCompose((PetscObject) dm, "dmAux", NULL);CHKERRQ(ierr);    ierr = PetscObjectCompose((PetscObject) dm, "A", NULL);CHKERRQ(ierr);  }  PetscFunctionReturn(0);}

/*@C    DMCompositeAddDM - adds a DM  vector to a DMComposite    Collective on DMComposite    Input Parameter:+    dm - the packer object-    dm - the DM object, if the DM is a da you will need to caste it with a (DM)    Level: advanced.seealso DMDestroy(), DMCompositeGather(), DMCompositeAddDM(), DMCreateGlobalVector(),         DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),         DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()@*/PetscErrorCode  DMCompositeAddDM(DM dmc,DM dm){  PetscErrorCode         ierr;  PetscInt               n,nlocal;  struct DMCompositeLink *mine,*next;  Vec                    global,local;  DM_Composite           *com = (DM_Composite*)dmc->data;  PetscFunctionBegin;  PetscValidHeaderSpecific(dmc,DM_CLASSID,1);  PetscValidHeaderSpecific(dm,DM_CLASSID,2);  next = com->next;  if (com->setup) SETERRQ(((PetscObject)dmc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot add a DM once you have used the DMComposite");  /* create new link */  ierr = PetscNew(struct DMCompositeLink,&mine);CHKERRQ(ierr);  ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);  ierr = DMGetGlobalVector(dm,&global);CHKERRQ(ierr);  ierr = VecGetLocalSize(global,&n);CHKERRQ(ierr);  ierr = DMRestoreGlobalVector(dm,&global);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm,&local);CHKERRQ(ierr);  ierr = VecGetSize(local,&nlocal);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm,&local);CHKERRQ(ierr);  mine->n      = n;  mine->nlocal = nlocal;  mine->dm     = dm;  mine->next   = PETSC_NULL;  com->n       += n;  /* add to end of list */  if (!next) {    com->next = mine;  } else {    while (next->next) next = next->next;    next->next = mine;  }  com->nDM++;  com->nmine++;  PetscFunctionReturn(0);}

static PetscErrorCode TestQ2Q1DA( void ){  DM             Q2_da,Q1_da,cda;  PetscInt       mx,my,mz;  Vec            coords,gcoords,gcoords2;  PetscErrorCode ierr;  mx=7;  my=11;  mz=13;  ierr=DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,mx,my,mz,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,3,2,0,0,0,&Q2_da);CHKERRQ(ierr);  ierr = DMDASetUniformCoordinates(Q2_da,-1.0,1.0,-2.0,2.0,-3.0,3.0);CHKERRQ(ierr);  ierr = DMDAGetCoordinates(Q2_da,&coords);CHKERRQ(ierr);  ierr = DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,mx,my,mz,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,3,1,0,0,0,&Q1_da);CHKERRQ(ierr);  ierr = DMDASetCoordinates(Q1_da,coords);CHKERRQ(ierr);  /* Get ghost coordinates one way */  ierr = DMDAGetGhostedCoordinates(Q1_da,&gcoords);CHKERRQ(ierr);  /* And another */  ierr = DMDAGetCoordinates(Q1_da,&coords);CHKERRQ(ierr);  ierr = DMDAGetCoordinateDA(Q1_da,&cda);CHKERRQ(ierr);  ierr = DMGetLocalVector(cda,&gcoords2);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(cda,coords,INSERT_VALUES,gcoords2);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(cda,coords,INSERT_VALUES,gcoords2);CHKERRQ(ierr);  ierr = CompareGhostedCoords(gcoords,gcoords2);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(cda,&gcoords2);CHKERRQ(ierr);  ierr = VecScale(coords,10.0);CHKERRQ(ierr);  ierr = VecScale(gcoords,10.0);CHKERRQ(ierr);  ierr = DMDAGetGhostedCoordinates(Q1_da,&gcoords2);CHKERRQ(ierr);  ierr = CompareGhostedCoords(gcoords,gcoords2);CHKERRQ(ierr);  ierr = DMDestroy(&Q2_da);CHKERRQ(ierr);  ierr = DMDestroy(&Q1_da);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode WaterSetInitialGuess(DM networkdm,Vec X){  PetscErrorCode ierr;  PetscInt       nv,ne;  const PetscInt *vtx,*edges;  Vec            localX;  PetscFunctionBegin;  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);  /* Get subnetwork info */  ierr = DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);CHKERRQ(ierr);  ierr = SetInitialGuess_Water(networkdm,localX,nv,ne,vtx,edges,NULL);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 FormIFunction(TS ts,PetscReal ptime,Vec X,Vec Xdot,Vec F,void *user){  DMDALocalInfo  info;  Field          **u,**udot,**fu;  PetscErrorCode ierr;  Vec            localX;  DM             da;  PetscFunctionBeginUser;  ierr = TSGetDM(ts,(DM*)&da);CHKERRQ(ierr);  ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr);  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 = DMDAVecGetArrayRead(da,localX,&u);CHKERRQ(ierr);  ierr = DMDAVecGetArrayRead(da,Xdot,&udot);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,F,&fu);CHKERRQ(ierr);  ierr = FormIFunctionLocal(&info,ptime,u,udot,fu,(AppCtx*)user);CHKERRQ(ierr);  ierr = DMDAVecRestoreArrayRead(da,localX,&u);CHKERRQ(ierr);  ierr = DMDAVecRestoreArrayRead(da,Xdot,&udot);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,F,&fu);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da,&localX);CHKERRQ(ierr);  PetscFunctionReturn(0);}