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

/*@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);}

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);}

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);}

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 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);}

/*  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 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);}

/*      Evaluates /integral_{-1}^{1} f*v_i  where v_i is the ith basis polynomial via the GLL nodes and weights, since the v_i      basis function is zero at all nodes except the ith one the integral is simply the weight_i * f(node_i)*/PetscErrorCode ComputeRhs(DM da,PetscGLL *gll,Vec b){  PetscErrorCode ierr;  PetscInt       i,j,xs,xn,n = gll->n;  PetscScalar    *bb,*xx;  PetscReal      xd;  Vec            blocal,xlocal;  PetscFunctionBegin;  ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xn,NULL,NULL);CHKERRQ(ierr);  xs   = xs/(n-1);  xn   = xn/(n-1);  ierr = DMGetLocalVector(da,&blocal);CHKERRQ(ierr);  ierr = VecZeroEntries(blocal);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,blocal,&bb);CHKERRQ(ierr);  ierr = DMGetCoordinatesLocal(da,&xlocal);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,xlocal,&xx);CHKERRQ(ierr);  /* loop over local spectral elements */  for (j=xs; j<xs+xn; j++) {    /* loop over GLL points in each element */    for (i=0; i<n; i++) {      xd              = xx[j*(n-1) + i];      bb[j*(n-1) + i] += -gll->weights[i]*(-20.*PETSC_PI*xd*PetscSinReal(5.*PETSC_PI*xd) + (2. - (5.*PETSC_PI)*(5.*PETSC_PI)*(xd*xd - 1.))*PetscCosReal(5.*PETSC_PI*xd));    }  }  ierr = DMDAVecRestoreArray(da,xlocal,&xx);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,blocal,&bb);CHKERRQ(ierr);  ierr = VecZeroEntries(b);CHKERRQ(ierr);  ierr = DMLocalToGlobalBegin(da,blocal,ADD_VALUES,b);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(da,blocal,ADD_VALUES,b);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da,&blocal);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);}

/*  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);}

/* 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);}

/*   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);}

/*  We integrate over each cell  (i, j+1)----(i+1, j+1)      | /         |      |  /        |      |   /       |      |    /      |      |     /     |      |      /    |      |       /   |  (i,   j)----(i+1, j)*/PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx){  UserContext    *user = (UserContext*)ctx;  PetscScalar    phi   = user->phi;  PetscScalar    *array;  PetscInt       ne,nc,i;  const PetscInt *e;  PetscErrorCode ierr;  Vec            blocal;  DM             da;  PetscFunctionBeginUser;  ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);  /* access a local vector with room for the ghost points */  ierr = DMGetLocalVector(da,&blocal);CHKERRQ(ierr);  ierr = VecGetArray(blocal, (PetscScalar**) &array);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 set each nodal value to phi (will actually do integration over element */    array[e[3*i]]   = phi;    array[e[3*i+1]] = phi;    array[e[3*i+2]] = phi;  }  ierr = VecRestoreArray(blocal, (PetscScalar**) &array);CHKERRQ(ierr);  ierr = DMDARestoreElements(da,&ne,&nc,&e);CHKERRQ(ierr);  /* add our partial sums over all processors into b */  ierr = DMLocalToGlobalBegin(da,blocal,ADD_VALUES,b);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(da,blocal, ADD_VALUES,b);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da,&blocal);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);}

/**   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);}

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 TestFieldProjection(DM dm, DM auxdm, DMLabel label, Vec la, const char name[], AppCtx *user){  PetscErrorCode (**afuncs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *);  void           (**funcs)(PetscInt, PetscInt, PetscInt,                           const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],                           const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],                           PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);  Vec               lx, lu;  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 = PetscMalloc2(Nf, &funcs, Nf, &afuncs);CHKERRQ(ierr);  for (f = 0; f < Nf; ++f) afuncs[f]  = linear;  funcs[0] = linear_vector;  funcs[1] = linear_scalar;  ierr = DMGetLocalVector(dm, &lu);CHKERRQ(ierr);  ierr = PetscStrcpy(lname, "Local Field Input ");CHKERRQ(ierr);  ierr = PetscStrcat(lname, name);CHKERRQ(ierr);  ierr = PetscObjectSetName((PetscObject) lu, lname);CHKERRQ(ierr);  if (!label) {ierr = DMProjectFunctionLocal(dm, 0.0, afuncs, NULL, INSERT_VALUES, lu);CHKERRQ(ierr);}  else        {ierr = DMProjectFunctionLabelLocal(dm, 0.0, label, 1, val, 0, NULL, afuncs, NULL, INSERT_VALUES, lu);CHKERRQ(ierr);}  ierr = VecViewFromOptions(lu, NULL, "-local_input_view");CHKERRQ(ierr);  ierr = DMGetLocalVector(dm, &lx);CHKERRQ(ierr);  ierr = PetscStrcpy(lname, "Local Field ");CHKERRQ(ierr);  ierr = PetscStrcat(lname, name);CHKERRQ(ierr);  ierr = PetscObjectSetName((PetscObject) lx, lname);CHKERRQ(ierr);  if (!label) {ierr = DMProjectFieldLocal(dm, 0.0, lu, funcs, INSERT_VALUES, lx);CHKERRQ(ierr);}  else        {ierr = DMProjectFieldLabelLocal(dm, 0.0, label, 1, val, 0, NULL, lu, funcs, INSERT_VALUES, lx);CHKERRQ(ierr);}  ierr = VecViewFromOptions(lx, NULL, "-local_field_view");CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm, &lx);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm, &lu);CHKERRQ(ierr);  ierr = PetscFree2(funcs, afuncs);CHKERRQ(ierr);  if (auxdm) {    ierr = PetscObjectCompose((PetscObject) dm, "dmAux", NULL);CHKERRQ(ierr);    ierr = PetscObjectCompose((PetscObject) dm, "A", NULL);CHKERRQ(ierr);  }  PetscFunctionReturn(0);}

PetscErrorCode OpForcing(Op op,DM dm,Vec F) {  PetscErrorCode ierr;  Vec X,Floc;  DM dmx;  const PetscScalar *x;  PetscScalar *f;  const PetscReal *B,*D,*w3;  PetscReal L[3];  PetscInt nelem,ne = op->ne,P,Q,P3,Q3;  PetscFunctionBegin;  ierr = PetscLogEventBegin(OP_Forcing,dm,F,0,0);CHKERRQ(ierr);  ierr = DMFEGetTensorEval(dm,&P,&Q,&B,&D,NULL,NULL,&w3);CHKERRQ(ierr);  P3 = P*P*P;  Q3 = Q*Q*Q;  ierr = DMFEGetUniformCoordinates(dm,L);CHKERRQ(ierr);  ierr = DMGetLocalVector(dm,&Floc);CHKERRQ(ierr);  ierr = DMGetCoordinateDM(dm,&dmx);CHKERRQ(ierr);  ierr = DMGetCoordinatesLocal(dm,&X);CHKERRQ(ierr);  ierr = DMFEGetNumElements(dm,&nelem);CHKERRQ(ierr);  ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);  ierr = VecZeroEntries(Floc);CHKERRQ(ierr);  ierr = VecGetArray(Floc,&f);CHKERRQ(ierr);  for (PetscInt e=0; e<nelem; e+=ne) {    PetscScalar fe[op->dof*P3*ne]_align,fq[op->dof][Q3][ne]_align,xe[3*P3*ne]_align,xq[3][Q3][ne]_align,dx[3][3][Q3][ne]_align,wdxdet[Q3][ne]_align;    ierr = DMFEExtractElements(dmx,x,e,ne,xe);CHKERRQ(ierr);    ierr = PetscMemzero(xq,sizeof xq);CHKERRQ(ierr);    ierr = TensorContract(op->Tensor3,B,B,B,TENSOR_EVAL,xe,xq[0][0]);CHKERRQ(ierr);    ierr = PetscMemzero(dx,sizeof dx);CHKERRQ(ierr);    ierr = TensorContract(op->Tensor3,D,B,B,TENSOR_EVAL,xe,dx[0][0][0]);CHKERRQ(ierr);    ierr = TensorContract(op->Tensor3,B,D,B,TENSOR_EVAL,xe,dx[1][0][0]);CHKERRQ(ierr);    ierr = TensorContract(op->Tensor3,B,B,D,TENSOR_EVAL,xe,dx[2][0][0]);CHKERRQ(ierr);    ierr = PointwiseJacobianInvert(ne,Q*Q*Q,w3,dx,wdxdet);CHKERRQ(ierr);    for (PetscInt i=0; i<Q3; i++) {      for (PetscInt l=0; l<ne; l++) {        PetscReal xx[] = {xq[0][i][l],xq[1][i][l],xq[2][i][l]};        PetscScalar fql[op->dof];        ierr = (op->PointwiseForcing)(op,xx,L,fql);CHKERRQ(ierr);        for (PetscInt d=0; d<op->dof; d++) fq[d][i][l] = wdxdet[i][l] * fql[d];      }    }    ierr = PetscMemzero(fe,sizeof fe);CHKERRQ(ierr);    ierr = TensorContract(op->TensorDOF,B,B,B,TENSOR_TRANSPOSE,fq[0][0],fe);CHKERRQ(ierr);    ierr = DMFESetElements(dm,f,e,ne,ADD_VALUES,DOMAIN_INTERIOR,fe);CHKERRQ(ierr);  }  ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);  ierr = VecRestoreArray(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);  ierr = PetscLogEventEnd(OP_Forcing,dm,F,0,0);CHKERRQ(ierr);  PetscFunctionReturn(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 FormFunctionMatlab(SNES snes,Vec X,Vec F,void *ptr){  AppCtx         *user = (AppCtx*)ptr;  PetscErrorCode ierr;  PetscInt       Mx,My;  PetscReal      lambda,hx,hy;  Vec            localX,localF;  MPI_Comm       comm;  DM             da;  PetscFunctionBeginUser;  ierr = SNESGetDM(snes,&da);CHKERRQ(ierr);  ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr);  ierr = DMGetLocalVector(da,&localF);CHKERRQ(ierr);  ierr = PetscObjectSetName((PetscObject)localX,"localX");CHKERRQ(ierr);  ierr = PetscObjectSetName((PetscObject)localF,"localF");CHKERRQ(ierr);  ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr);  lambda = user->param;  hx     = 1.0/(PetscReal)(Mx-1);  hy     = 1.0/(PetscReal)(My-1);  ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr);  /*     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);  ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localX);CHKERRQ(ierr);  ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(comm),"localF=ex5m(localX,%18.16e,%18.16e,%18.16e)",hx,hy,lambda);CHKERRQ(ierr);  ierr = PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localF);CHKERRQ(ierr);  /*     Insert values into global vector  */  ierr = DMLocalToGlobalBegin(da,localF,INSERT_VALUES,F);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(da,localF,INSERT_VALUES,F);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da,&localX);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(da,&localF);CHKERRQ(ierr);  PetscFunctionReturn(0);}

/* no discretization is given so DMGetNumFields yields 0 */static PetscErrorCode test0(DM dm, AppCtx *options){  Vec            locX;  PetscErrorCode ierr;  PetscFunctionBegin;  ierr = DMGetLocalVector(dm, &locX);CHKERRQ(ierr);  ierr = DMRestoreLocalVector(dm, &locX);CHKERRQ(ierr);  PetscFunctionReturn(0);}

int main(int argc, char **argv){  DM             dm, subdm, auxdm;  Vec            la;  AppCtx         user;  PetscErrorCode ierr;  ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;  ierr = ProcessOptions(&user);CHKERRQ(ierr);  ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr);  ierr = SetupDiscretization(dm, user.dim, user.cellSimplex, &user);CHKERRQ(ierr);  /* Volumetric Mesh Projection */  ierr = TestFunctionProjection(dm, NULL, NULL, NULL, "Volumetric Primary", &user);CHKERRQ(ierr);  ierr = TestFieldProjection(dm, NULL, NULL, NULL, "Volumetric Primary", &user);CHKERRQ(ierr);  if (user.auxfield) {    /* Volumetric Mesh Projection with Volumetric Data */    ierr = CreateAuxiliaryData(dm, &auxdm, &la, &user);CHKERRQ(ierr);    ierr = TestFunctionProjection(dm, auxdm, NULL, la, "Volumetric Primary and Volumetric Auxiliary", &user);CHKERRQ(ierr);    ierr = TestFieldProjection(dm, auxdm, NULL, la, "Volumetric Primary and Volumetric Auxiliary", &user);CHKERRQ(ierr);    ierr = VecDestroy(&la);CHKERRQ(ierr);    /* Update of Volumetric Auxiliary Data with primary Volumetric Data */    ierr = DMGetLocalVector(dm, &la);CHKERRQ(ierr);    ierr = VecSet(la, 1.0);CHKERRQ(ierr);    ierr = TestFieldProjection(auxdm, dm, NULL, la, "Volumetric Auxiliary Update with Volumetric Primary", &user);CHKERRQ(ierr);    ierr = DMRestoreLocalVector(dm, &la);CHKERRQ(ierr);    ierr = DMDestroy(&auxdm);CHKERRQ(ierr);  }  if (user.submesh) {    DMLabel bdLabel;    /* Boundary Mesh Projection */    ierr = CreateBoundaryMesh(dm, &bdLabel, &subdm, &user);CHKERRQ(ierr);    ierr = TestFunctionProjection(subdm, NULL, NULL, NULL, "Boundary Primary", &user);CHKERRQ(ierr);    ierr = TestFieldProjection(subdm, NULL, NULL, NULL, "Boundary Primary", &user);CHKERRQ(ierr);    if (user.auxfield) {      /* Boundary Mesh Projection with Boundary Data */      ierr = CreateAuxiliaryData(subdm, &auxdm, &la, &user);CHKERRQ(ierr);      ierr = TestFunctionProjection(subdm, auxdm, NULL, la, "Boundary Primary and Boundary Auxiliary", &user);CHKERRQ(ierr);      ierr = TestFieldProjection(subdm, auxdm, NULL, la, "Boundary Primary and Boundary Auxiliary", &user);CHKERRQ(ierr);      ierr = VecDestroy(&la);CHKERRQ(ierr);      ierr = DMDestroy(&auxdm);CHKERRQ(ierr);      /* Volumetric Mesh Projection with Boundary Data */      ierr = CreateAuxiliaryData(subdm, &auxdm, &la, &user);CHKERRQ(ierr);      ierr = TestFunctionProjection(dm, auxdm, bdLabel, la, "Volumetric Primary and Boundary Auxiliary", &user);CHKERRQ(ierr);      ierr = TestFieldProjection(dm, auxdm, bdLabel, la, "Volumetric Primary and Boundary Auxiliary", &user);CHKERRQ(ierr);      ierr = VecDestroy(&la);CHKERRQ(ierr);      ierr = DMDestroy(&auxdm);CHKERRQ(ierr);    }    ierr = DMLabelDestroy(&bdLabel);CHKERRQ(ierr);    ierr = DMDestroy(&subdm);CHKERRQ(ierr);  }  ierr = DMDestroy(&dm);CHKERRQ(ierr);  ierr = PetscFinalize();  return ierr;}

/* no discretization is given so DMGetNumFields and PetscDSGetTotalDimension yield 0 */static PetscErrorCode test1(DM dm, AppCtx *options){  IS             cells;  Vec            locX, locX_t, locA;  PetscScalar    *u, *u_t, *a;  PetscErrorCode ierr;  PetscFunctionBegin;  ierr = ISCreateStride(PETSC_COMM_SELF, 0, 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);}

/*     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 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);}

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);}