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

int main(int argc,char **argv){  PetscInt       M = 14,dof = 1,s = 1,ratio = 2,dim = 2;  PetscErrorCode ierr;  DM             da_c,da_f;  Vec            v_c,v_f;  Mat            I;  PetscScalar    one = 1.0;  MPI_Comm       comm_f, comm_c;  ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);   ierr = PetscOptionsGetInt(PETSC_NULL,"-dim",&dim,PETSC_NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(PETSC_NULL,"-sw",&s,PETSC_NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(PETSC_NULL,"-ratio",&ratio,PETSC_NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(PETSC_NULL,"-dof",&dof,PETSC_NULL);CHKERRQ(ierr);  comm_f = PETSC_COMM_WORLD;  ierr = DMDASplitComm2d(comm_f,M,M,s,&comm_c);CHKERRQ(ierr);      /* Set up the array */   if (dim == 2) {    ierr = DMDACreate2d(comm_c,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,M,M,PETSC_DECIDE,PETSC_DECIDE,dof,s,PETSC_NULL,PETSC_NULL,&da_c);CHKERRQ(ierr);    M    = ratio*(M-1) + 1;    ierr = DMDACreate2d(comm_f,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,M,M,PETSC_DECIDE,PETSC_DECIDE,dof,s,PETSC_NULL,PETSC_NULL,&da_f);CHKERRQ(ierr);  } else if (dim == 3) {    ;  }  ierr = DMCreateGlobalVector(da_c,&v_c);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da_f,&v_f);CHKERRQ(ierr);  ierr = VecSet(v_c,one);CHKERRQ(ierr);  ierr = DMCreateInterpolation(da_c,da_f,&I,PETSC_NULL);CHKERRQ(ierr);  ierr = MatInterpolate(I,v_c,v_f);CHKERRQ(ierr);  ierr = VecView(v_f,PETSC_VIEWER_STDOUT_(comm_f));CHKERRQ(ierr);  ierr = MatRestrict(I,v_f,v_c);CHKERRQ(ierr);  ierr = VecView(v_c,PETSC_VIEWER_STDOUT_(comm_c));CHKERRQ(ierr);  ierr = MatDestroy(&I);CHKERRQ(ierr);  ierr = VecDestroy(&v_c);CHKERRQ(ierr);  ierr = DMDestroy(&da_c);CHKERRQ(ierr);  ierr = VecDestroy(&v_f);CHKERRQ(ierr);  ierr = DMDestroy(&da_f);CHKERRQ(ierr);  ierr = PetscFinalize();  return 0;}

int main(int argc,char **argv){  PetscInt       dof = 2,M = 3,N = 3,m = PETSC_DECIDE,n = PETSC_DECIDE;  PetscErrorCode ierr;  DM             da;  Vec            global,local;  ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(0,"-dof",&dof,0);CHKERRQ(ierr);  /* Create distributed array and get vectors */  ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,                    M,N,m,n,dof,1,PETSC_NULL,PETSC_NULL,&da);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da,&global);CHKERRQ(ierr);  ierr = DMCreateLocalVector(da,&local);CHKERRQ(ierr);  ierr = doit(da,global);CHKERRQ(ierr);  ierr = VecView(global,0);CHKERRQ(ierr);  /* Free memory */  ierr = VecDestroy(&local);CHKERRQ(ierr);  ierr = VecDestroy(&global);CHKERRQ(ierr);  ierr = DMDestroy(&da);CHKERRQ(ierr);  ierr = PetscFinalize();  return 0;}

inline voidDMCreate<2         >(MPI_Comm const&                       comm,           DMBoundaryTypeVector<2> const&        boundaryTypes,           DMDAStencilType const&                stencilType,           VectorDi<2> const&                    globalSize,           VectorDi<2> const&                    processorSize,           int const&                            dof,           int const&                            stencilWidth,           VectorDConstPetscIntPointer<2> const& localSizes,           DM*                                   da) {  DMDACreate2d(comm,               boundaryTypes(0),               boundaryTypes(1),               stencilType,               globalSize(0),               globalSize(1),               processorSize(0),               processorSize(1),               dof,               stencilWidth,               localSizes(0).get(),               localSizes(1).get(),               da);}

int main(int argc,char **argv){  KSP            ksp;  DM             da;  UserContext    user;  PetscInt       bc;  PetscErrorCode ierr;  PetscInitialize(&argc,&argv,(char *)0,help);  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);  ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-11,-11,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,&da);CHKERRQ(ierr);  ierr = KSPSetDM(ksp,(DM)da);  ierr = DMSetApplicationContext(da,&user);CHKERRQ(ierr);  user.uu = 1.0;  user.tt = 1.0;  bc   = (PetscInt)NEUMANN; // Use Neumann Boundary Conditions  user.bcType = (BCType)bc;  ierr = KSPSetComputeRHS(ksp,ComputeRHS,&user);CHKERRQ(ierr);  ierr = KSPSetComputeOperators(ksp,ComputeJacobian,&user);CHKERRQ(ierr);  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);  ierr = KSPSolve(ksp,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);  ierr = DMDestroy(&da);CHKERRQ(ierr);  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);  ierr = PetscFinalize();CHKERRQ(ierr);  return 0;}

int main(int argc,char **argv){  SNES           snes;  PetscErrorCode ierr;  PetscInt       its,lits;  PetscReal      litspit;  DM             da;  PetscInitialize(&argc,&argv,PETSC_NULL,help);  /*      Set the DMDA (grid structure) for the grids.  */  ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-5,-5,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);CHKERRQ(ierr);  ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,PETSC_NULL);CHKERRQ(ierr);  ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);  ierr = SNESSetDM(snes,da);CHKERRQ(ierr);  ierr = DMDestroy(&da);CHKERRQ(ierr);  ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);  ierr = SNESSolve(snes,0,0);CHKERRQ(ierr);  ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr);  ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr);  litspit = ((PetscReal)lits)/((PetscReal)its);  ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D/n",its);CHKERRQ(ierr);  ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D/n",lits);CHKERRQ(ierr);  ierr = PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / SNES = %e/n",litspit);CHKERRQ(ierr);  ierr = SNESDestroy(&snes);CHKERRQ(ierr);  ierr = PetscFinalize();  return 0;}

int main(int argc,char **argv){    KSP            ksp;    DM             da;    UserContext    user;    const char     *bcTypes[2] = {"dirichlet","neumann"};    PetscErrorCode ierr;    PetscInt       bc;    Vec            b,x;    PetscInitialize(&argc,&argv,(char *)0,help);    ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);    CHKERRQ(ierr);    ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-3,-3,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);    CHKERRQ(ierr);    ierr = DMDASetUniformCoordinates(da,0,1,0,1,0,0);    CHKERRQ(ierr);    ierr = DMDASetFieldName(da,0,"Pressure");    CHKERRQ(ierr);    ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "DMqq");    user.rho    = 1.0;    ierr        = PetscOptionsReal("-rho", "The conductivity", "ex29.c", user.rho, &user.rho, PETSC_NULL);    CHKERRQ(ierr);    user.nu     = 0.1;    ierr        = PetscOptionsReal("-nu", "The width of the Gaussian source", "ex29.c", user.nu, &user.nu, PETSC_NULL);    CHKERRQ(ierr);    bc          = (PetscInt)DIRICHLET;    ierr        = PetscOptionsEList("-bc_type","Type of boundary condition","ex29.c",bcTypes,2,bcTypes[0],&bc,PETSC_NULL);    CHKERRQ(ierr);    user.bcType = (BCType)bc;    ierr = PetscOptionsEnd();    ierr = KSPSetComputeRHS(ksp,ComputeRHS,&user);    CHKERRQ(ierr);    ierr = KSPSetComputeOperators(ksp,ComputeMatrix,&user);    CHKERRQ(ierr);    ierr = KSPSetDM(ksp,da);    CHKERRQ(ierr);    ierr = KSPSetFromOptions(ksp);    CHKERRQ(ierr);    ierr = KSPSetUp(ksp);    CHKERRQ(ierr);    ierr = KSPSolve(ksp,PETSC_NULL,PETSC_NULL);    CHKERRQ(ierr);    ierr = KSPGetSolution(ksp,&x);    CHKERRQ(ierr);    ierr = KSPGetRhs(ksp,&b);    CHKERRQ(ierr);    ierr = DMDestroy(&da);    CHKERRQ(ierr);    ierr = KSPDestroy(&ksp);    CHKERRQ(ierr);    ierr = PetscFinalize();    return 0;}

PETSC_EXTERN void PETSC_STDCALL dmdacreate2d_(MPI_Comm *comm,DMBoundaryType *bx,DMBoundaryType *by,DMDAStencilType                  *stencil_type,PetscInt *M,PetscInt *N,PetscInt *m,PetscInt *n,PetscInt *w,                  PetscInt *s,PetscInt *lx,PetscInt *ly,DM *inra,PetscErrorCode *ierr){  CHKFORTRANNULLINTEGER(lx);  CHKFORTRANNULLINTEGER(ly);  *ierr = DMDACreate2d(MPI_Comm_f2c(*(MPI_Fint*)&*comm),*bx,*by,*stencil_type,*M,*N,*m,*n,*w,*s,lx,ly,inra);}

int main(int argc,char **argv){  PetscErrorCode ierr;  DM             da2D;  PetscInt       i,j,ixs, ixm, iys, iym;;  PetscViewer    H5viewer;  PetscScalar    xm=-1.0, xp=1.0;  PetscScalar    ym=-1.0, yp=1.0;  PetscScalar    value=1.0,dx,dy;  PetscInt       Nx=40, Ny=40;  Vec            gauss;  PetscScalar    **gauss_ptr;  dx=(xp-xm)/(Nx-1);  dy=(yp-ym)/(Ny-1);  // Initialize the Petsc context  ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);  // Build of the DMDA  ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,Nx,Ny,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,&da2D);CHKERRQ(ierr);  // Set the coordinates  DMDASetUniformCoordinates(da2D, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);  // Declare gauss as a DMDA component  ierr = DMCreateGlobalVector(da2D,&gauss);CHKERRQ(ierr);  ierr = PetscObjectSetName((PetscObject) gauss, "pressure");CHKERRQ(ierr);  // Initialize vector gauss with a constant value (=1)  ierr = VecSet(gauss,value);CHKERRQ(ierr);  // Get the coordinates of the corners for each process  ierr = DMDAGetCorners(da2D, &ixs, &iys, 0, &ixm, &iym, 0);CHKERRQ(ierr);  /* Build the gaussian profile (exp(-x^2-y^2)) */  ierr = DMDAVecGetArray(da2D,gauss,&gauss_ptr);CHKERRQ(ierr);  for (j=iys; j<iys+iym; j++){    for (i=ixs; i<ixs+ixm; i++){      gauss_ptr[j][i]=exp(-(xm+i*dx)*(xm+i*dx)-(ym+j*dy)*(ym+j*dy));    }  }  ierr = DMDAVecRestoreArray(da2D,gauss,&gauss_ptr);CHKERRQ(ierr);  // Create the HDF5 viewer  ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD,"gauss.h5",FILE_MODE_WRITE,&H5viewer);CHKERRQ(ierr);  // Write the H5 file  ierr = VecView(gauss,H5viewer);CHKERRQ(ierr);  // Cleaning stage  ierr = PetscViewerDestroy(&H5viewer);CHKERRQ(ierr);  ierr = VecDestroy(&gauss);CHKERRQ(ierr);  ierr = DMDestroy(&da2D);CHKERRQ(ierr);  ierr = PetscFinalize();    return 0;}

PetscErrorCode createDomain(DM *dm, int nx, int ny){  PetscErrorCode ierr;  PetscFunctionBegin;  ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,                       DMDA_STENCIL_STAR,                      nx, ny, PETSC_DECIDE, PETSC_DECIDE,                      1, 1, PETSC_NULL, PETSC_NULL, dm);CHKERRQ(ierr);  PetscFunctionReturn(0);}

void createArrays(simInfo *data){	PetscErrorCode   ierr;	PetscInt         *lx, *ly, m, n;	const PetscInt   *lxu, *lyu;		// Create distributed array and get vectors	ierr = PetscPrintf(PETSC_COMM_WORLD, "/nCreate staggered U, V and P vectors [DMDACreate2d, DMDAGetInfo, DMDAGetOwnershipRanges, DMCreateGlobalVector]"); CHKERRV(ierr);    hline();	ierr = DMCompositeCreate(PETSC_COMM_WORLD, &(data->pack)); CHKERRV(ierr);		// create u DA	ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_GHOSTED, DMDA_STENCIL_BOX, data->nx-1, data->ny, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &(data->uda)); CHKERRV(ierr);	ierr = DMCompositeAddDM(data->pack, data->uda); CHKERRV(ierr);		// determine process distribution of v	ierr = DMDAGetInfo(data->uda, NULL, NULL, NULL, NULL, &m, &n, NULL, NULL, NULL, NULL, NULL, NULL, NULL); CHKERRV(ierr);	ierr = DMDAGetOwnershipRanges(data->uda, &lxu, &lyu, NULL); CHKERRV(ierr);	ierr = PetscMalloc(m*sizeof(*lx), &lx); CHKERRV(ierr);	ierr = PetscMalloc(n*sizeof(*ly), &ly); CHKERRV(ierr);	ierr = PetscMemcpy(lx ,lxu, m*sizeof(*lx)); CHKERRV(ierr);	ierr = PetscMemcpy(ly ,lyu, n*sizeof(*ly)); CHKERRV(ierr);	lx[m-1]++;		// create v DA	ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_GHOSTED, DMDA_STENCIL_STAR, data->nx, data->ny, m, n, 1, 1, lx, ly, &(data->pda)); CHKERRV(ierr);		// create vectors p and bc2	ierr = DMCreateGlobalVector(data->pda, &(data->pGlobal)); CHKERRV(ierr);	ierr = VecDuplicate(data->pGlobal, &(data->bc2)); CHKERRV(ierr);	ly[n-1]--;		// create p DA	ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_GHOSTED, DMDA_STENCIL_BOX, data->nx, data->ny-1, m, n, 1, 1, lx, ly, &(data->vda)); CHKERRV(ierr);	ierr = DMCompositeAddDM(data->pack, data->vda); CHKERRV(ierr);		PetscFree(lx);	PetscFree(ly);		// create vector uPacked	ierr = DMCreateGlobalVector(data->pack, &(data->uPacked)); CHKERRV(ierr);}

PetscErrorCodePETScExternalSolverCreate(MPI_Comm comm, TS * ts){  DM da;  PetscErrorCode ierr;  PetscFunctionBeginUser;  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Create distributed array (DMDA) to manage parallel grid and vectors  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = DMDACreate2d(comm,                      DM_BOUNDARY_NONE,                      DM_BOUNDARY_NONE,                      DMDA_STENCIL_STAR,                      11,                      11,                      PETSC_DECIDE,                      PETSC_DECIDE,                      1,                      1,                      NULL,                      NULL,                      &da);  CHKERRQ(ierr);  ierr = DMSetFromOptions(da);  CHKERRQ(ierr);  ierr = DMSetUp(da);  CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -     Create timestepping solver context     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = TSCreate(comm, ts);  CHKERRQ(ierr);  ierr = TSSetProblemType(*ts, TS_NONLINEAR);  CHKERRQ(ierr);  ierr = TSSetType(*ts, TSBEULER);  CHKERRQ(ierr);  ierr = TSSetDM(*ts, da);  CHKERRQ(ierr);  ierr = DMDestroy(&da);  CHKERRQ(ierr);  ierr = TSSetIFunction(*ts, NULL, FormIFunction, nullptr);  CHKERRQ(ierr);  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -   Set Jacobian evaluation routine  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  ierr = TSSetIJacobian(*ts, NULL, NULL, FormIJacobian, NULL);  CHKERRQ(ierr);  PetscFunctionReturn(0);}

int main(int argc, char **argv){  PetscErrorCode ierr;  Vec            x;  /* Solution vector */  TS             ts;   /* Time-stepping context */  AppCtx         user; /* Application context */  Mat            J;  PetscViewer    viewer;  PetscInitialize(&argc,&argv,"petscopt_ex6", help);  /* Get physics and time parameters */  ierr = Parameter_settings(&user);CHKERRQ(ierr);  /* Create a 2D DA with dof = 1 */  ierr = DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&user.da);CHKERRQ(ierr);  /* Set x and y coordinates */  ierr = DMDASetUniformCoordinates(user.da,user.xmin,user.xmax,user.ymin,user.ymax,NULL,NULL);CHKERRQ(ierr);  /* Get global vector x from DM  */  ierr = DMCreateGlobalVector(user.da,&x);CHKERRQ(ierr);  ierr = ini_bou(x,&user);CHKERRQ(ierr);  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"ini_x",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);  ierr = VecView(x,viewer);CHKERRQ(ierr);  ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);  /* Get Jacobian matrix structure from the da */  ierr = DMSetMatType(user.da,MATAIJ);CHKERRQ(ierr);  ierr = DMCreateMatrix(user.da,&J);CHKERRQ(ierr);  ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);  ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);  ierr = TSSetIFunction(ts,NULL,IFunction,&user);CHKERRQ(ierr);  ierr = TSSetIJacobian(ts,J,J,IJacobian,&user);CHKERRQ(ierr);  ierr = TSSetApplicationContext(ts,&user);CHKERRQ(ierr);  ierr = TSSetDuration(ts,PETSC_DEFAULT,user.tmax);CHKERRQ(ierr);  ierr = TSSetInitialTimeStep(ts,user.t0,.005);CHKERRQ(ierr);  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);  ierr = TSSetPostStep(ts,PostStep);CHKERRQ(ierr);  ierr = TSSolve(ts,x);CHKERRQ(ierr);  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"fin_x",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);  ierr = VecView(x,viewer);CHKERRQ(ierr);  ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);  ierr = VecDestroy(&x);CHKERRQ(ierr);  ierr = MatDestroy(&J);CHKERRQ(ierr);  ierr = DMDestroy(&user.da);CHKERRQ(ierr);  ierr = TSDestroy(&ts);CHKERRQ(ierr);  PetscFinalize();  return 0;}

int main(int argc, char *argv[]){  DM             da, daX, daY;  DMDALocalInfo  info;  MPI_Comm       commX, commY;  Vec            basisX, basisY;  PetscScalar    **arrayX, **arrayY;  const PetscInt *lx, *ly;  PetscInt       M     = 3, N = 3;  PetscInt       p     = 1;  PetscInt       numGP = 3;  PetscInt       dof   = 2*(p+1)*numGP;  PetscMPIInt    rank, subsize, subrank;  PetscErrorCode ierr;  ierr = PetscInitialize(&argc,&argv,0,help);CHKERRQ(ierr);  ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr);  /* Create 2D DMDA */  ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da);CHKERRQ(ierr);  /* Create 1D DMDAs along two directions */  ierr = DMDAGetOwnershipRanges(da, &lx, &ly, NULL);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(da, &info);CHKERRQ(ierr);  ierr = DMDAGetProcessorSubsets(da, DMDA_X, &commX);CHKERRQ(ierr);  ierr = DMDAGetProcessorSubsets(da, DMDA_Y, &commY);CHKERRQ(ierr);  ierr = MPI_Comm_size(commX, &subsize);CHKERRQ(ierr);  ierr = MPI_Comm_rank(commX, &subrank);CHKERRQ(ierr);  ierr = PetscPrintf(PETSC_COMM_SELF, "[%d]X subrank: %d subsize: %d/n", rank, subrank, subsize);  ierr = MPI_Comm_size(commY, &subsize);CHKERRQ(ierr);  ierr = MPI_Comm_rank(commY, &subrank);CHKERRQ(ierr);  ierr = PetscPrintf(PETSC_COMM_SELF, "[%d]Y subrank: %d subsize: %d/n", rank, subrank, subsize);  ierr = DMDACreate1d(commX, DM_BOUNDARY_NONE, M, dof, 1, lx, &daX);CHKERRQ(ierr);  ierr = DMDACreate1d(commY, DM_BOUNDARY_NONE, N, dof, 1, ly, &daY);CHKERRQ(ierr);  /* Create 1D vectors for basis functions */  ierr = DMGetGlobalVector(daX, &basisX);CHKERRQ(ierr);  ierr = DMGetGlobalVector(daY, &basisY);CHKERRQ(ierr);  /* Extract basis functions */  ierr = DMDAVecGetArrayDOF(daX, basisX, &arrayX);CHKERRQ(ierr);  ierr = DMDAVecGetArrayDOF(daY, basisY, &arrayY);CHKERRQ(ierr);  /*arrayX[i][ndof]; */  /*arrayY[j][ndof]; */  ierr = DMDAVecRestoreArrayDOF(daX, basisX, &arrayX);CHKERRQ(ierr);  ierr = DMDAVecRestoreArrayDOF(daY, basisY, &arrayY);CHKERRQ(ierr);  /* Return basis vectors */  ierr = DMRestoreGlobalVector(daX, &basisX);CHKERRQ(ierr);  ierr = DMRestoreGlobalVector(daY, &basisY);CHKERRQ(ierr);  /* Cleanup */  ierr = DMDestroy(&daX);CHKERRQ(ierr);  ierr = DMDestroy(&daY);CHKERRQ(ierr);  ierr = DMDestroy(&da);CHKERRQ(ierr);  ierr = PetscFinalize();  return 0;}

/*   Sets up a monitor that will display He as a function of space and cluster size for each time step*/PetscErrorCode MyMonitorSetUp(TS ts){  DM             da;  PetscErrorCode ierr;  PetscInt       xi,xs,xm,*idx,M,xj,cnt = 0,dof = 3*N + N*N;  const PetscInt *lx;  Vec            C;  MyMonitorCtx   *ctx;  PetscBool      flg;  IS             is;  char           ycoor[32];  PetscReal      valuebounds[4]  = {0, 1.2, 0, 1.2};  PetscFunctionBeginUser;  ierr = PetscOptionsHasName(PETSC_NULL,"-mymonitor",&flg);CHKERRQ(ierr);  if (!flg) PetscFunctionReturn(0);  ierr = TSGetDM(ts,&da);CHKERRQ(ierr);  ierr = PetscNew(MyMonitorCtx,&ctx);CHKERRQ(ierr);  ierr = PetscViewerDrawOpen(((PetscObject)da)->comm,PETSC_NULL,"",PETSC_DECIDE,PETSC_DECIDE,600,400,&ctx->viewer);CHKERRQ(ierr);  ierr = DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);  ierr = DMDAGetInfo(da,PETSC_IGNORE,&M,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr);  ierr = DMDAGetOwnershipRanges(da,&lx,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);  ierr = DMDACreate2d(((PetscObject)da)->comm,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,M,N,PETSC_DETERMINE,1,2,1,lx,PETSC_NULL,&ctx->da);CHKERRQ(ierr);  ierr = DMDASetFieldName(ctx->da,0,"He");CHKERRQ(ierr);  ierr = DMDASetFieldName(ctx->da,1,"V");CHKERRQ(ierr);  ierr = DMDASetCoordinateName(ctx->da,0,"X coordinate direction");CHKERRQ(ierr);  ierr = PetscSNPrintf(ycoor,32,"%D ... Cluster size ... 1",N);CHKERRQ(ierr);  ierr = DMDASetCoordinateName(ctx->da,1,ycoor);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(ctx->da,&ctx->He);CHKERRQ(ierr);  ierr = PetscMalloc(2*N*xm*sizeof(PetscInt),&idx);CHKERRQ(ierr);  cnt  = 0;  for (xj=0; xj<N; xj++) {    for (xi=xs; xi<xs+xm; xi++) {      idx[cnt++] = dof*xi + xj;      idx[cnt++] = dof*xi + xj + N;    }  }  ierr = ISCreateGeneral(((PetscObject)ts)->comm,2*N*xm,idx,PETSC_OWN_POINTER,&is);CHKERRQ(ierr);  ierr = TSGetSolution(ts,&C);CHKERRQ(ierr);  ierr = VecScatterCreate(C,is,ctx->He,PETSC_NULL,&ctx->scatter);CHKERRQ(ierr);  ierr = ISDestroy(&is);CHKERRQ(ierr);  /* sets the bounds on the contour plot values so the colors mean the same thing for different timesteps */  ierr = PetscViewerDrawSetBounds(ctx->viewer,2,valuebounds);CHKERRQ(ierr);  ierr = TSMonitorSet(ts,MyMonitorMonitor,ctx,MyMonitorDestroy);CHKERRQ(ierr);  PetscFunctionReturn(0);}

int main(int argc,char **argv){  PetscInt         M = 5,N = 4,P = 3, m = PETSC_DECIDE,n = PETSC_DECIDE,p = PETSC_DECIDE,dim = 1;  PetscErrorCode   ierr;  DM               dac,daf;  DMDABoundaryType bx=DMDA_BOUNDARY_NONE,by=DMDA_BOUNDARY_NONE,bz=DMDA_BOUNDARY_NONE;  DMDAStencilType  stype = DMDA_STENCIL_BOX;  Mat              A;  ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);  /* Read options */  ierr = PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(PETSC_NULL,"-N",&N,PETSC_NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(PETSC_NULL,"-P",&P,PETSC_NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(PETSC_NULL,"-p",&p,PETSC_NULL);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(PETSC_NULL,"-dim",&dim,PETSC_NULL);CHKERRQ(ierr);  /* Create distributed array and get vectors */  if (dim == 1) {    ierr = DMDACreate1d(PETSC_COMM_WORLD,bx,M,1,1,PETSC_NULL,&dac);CHKERRQ(ierr);  } else if (dim == 2) {    ierr = DMDACreate2d(PETSC_COMM_WORLD,bx,by,stype,M,N,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,&dac);CHKERRQ(ierr);  } else if (dim == 3) {    ierr = DMDACreate3d(PETSC_COMM_WORLD,bx,by,bz,stype,M,N,P,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,PETSC_NULL,&dac);CHKERRQ(ierr);  }  ierr = DMRefine(dac,PETSC_COMM_WORLD,&daf);CHKERRQ(ierr);  ierr = DMDASetUniformCoordinates(dac,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr);  if (dim == 1) {    ierr = SetCoordinates1d(daf);CHKERRQ(ierr);  } else if (dim == 2) {    ierr = SetCoordinates2d(daf);CHKERRQ(ierr);  } else if (dim == 3) {    ierr = SetCoordinates3d(daf);CHKERRQ(ierr);  }  ierr = DMCreateInterpolation(dac,daf,&A,0);CHKERRQ(ierr);  ierr = MatViewFromOptions(A,"-mat_view");CHKERRQ(ierr);  /* Free memory */  ierr = DMDestroy(&dac);CHKERRQ(ierr);  ierr = DMDestroy(&daf);CHKERRQ(ierr);  ierr = MatDestroy(&A);CHKERRQ(ierr);  ierr = PetscFinalize();  return 0;}

int main(int argc, char **argv){  PetscErrorCode ierr;  Vec            x;  /* Solution vector */  TS             ts;   /* Time-stepping context */  AppCtx         user; /* Application context */  PetscViewer    viewer;  PetscInitialize(&argc,&argv,"petscopt_ex7", help);  /* Get physics and time parameters */  ierr = Parameter_settings(&user);CHKERRQ(ierr);  /* Create a 2D DA with dof = 1 */  ierr = DMDACreate2d(PETSC_COMM_WORLD,user.bx,user.by,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,1,user.st_width,NULL,NULL,&user.da);CHKERRQ(ierr);  /* Set x and y coordinates */  ierr = DMDASetUniformCoordinates(user.da,user.xmin,user.xmax,user.ymin,user.ymax,0,0);CHKERRQ(ierr);  ierr = DMDASetCoordinateName(user.da,0,"X - the angle");  ierr = DMDASetCoordinateName(user.da,1,"Y - the speed");  /* Get global vector x from DM  */  ierr = DMCreateGlobalVector(user.da,&x);CHKERRQ(ierr);  ierr = ini_bou(x,&user);CHKERRQ(ierr);  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"ini_x",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);  ierr = VecView(x,viewer);CHKERRQ(ierr);  ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);  ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);  ierr = TSSetDM(ts,user.da);CHKERRQ(ierr);  ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);  ierr = TSSetType(ts,TSARKIMEX);CHKERRQ(ierr);  ierr = TSSetIFunction(ts,NULL,IFunction,&user);CHKERRQ(ierr);  /*  ierr = TSSetIJacobian(ts,NULL,NULL,IJacobian,&user);CHKERRQ(ierr);  */  ierr = TSSetApplicationContext(ts,&user);CHKERRQ(ierr);  ierr = TSSetInitialTimeStep(ts,0.0,.005);CHKERRQ(ierr);  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);  ierr = TSSetPostStep(ts,PostStep);CHKERRQ(ierr);  ierr = TSSolve(ts,x);CHKERRQ(ierr);  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"fin_x",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);  ierr = VecView(x,viewer);CHKERRQ(ierr);  ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);  ierr = VecDestroy(&x);CHKERRQ(ierr);  ierr = DMDestroy(&user.da);CHKERRQ(ierr);  ierr = TSDestroy(&ts);CHKERRQ(ierr);  PetscFinalize();  return 0;}

int main(int argc,char **argv){  PetscErrorCode ierr;  PetscInt       M = -8, N = -8,stencil_width = 1, dof = 1,m,n,xstart,ystart,i,j,c;  DM             da;  Vec            global,local;  PetscScalar    ***vglobal;  ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);  PetscFunctionBeginUser;  ierr = PetscOptionsGetInt(0,"-stencil_width",&stencil_width,0);CHKERRQ(ierr);  ierr = PetscOptionsGetInt(0,"-dof",&dof,0);CHKERRQ(ierr);  ierr = DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_MIRROR,DMDA_BOUNDARY_MIRROR,DMDA_STENCIL_STAR,M,N,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,&da);CHKERRQ(ierr);  ierr = DMDAGetCorners(da,&xstart,&ystart,0,&m,&n,0);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da,&global);CHKERRQ(ierr);  ierr = DMDAVecGetArrayDOF(da,global,&vglobal);CHKERRQ(ierr);  for (j=ystart; j<ystart+n; j++) {    for (i=xstart; i<xstart+m; i++) {      for (c=0; c<dof; c++) {        vglobal[j][i][c] = 100*j + 10*(i+1) + c;      }    }  }  ierr = DMDAVecRestoreArrayDOF(da,global,&vglobal);CHKERRQ(ierr);  ierr = DMCreateLocalVector(da,&local);CHKERRQ(ierr);  ierr = DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);CHKERRQ(ierr);  ierr = DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);CHKERRQ(ierr);  ierr = PetscSequentialPhaseBegin(PETSC_COMM_WORLD,1);CHKERRQ(ierr);  ierr = VecView(local,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);  ierr = PetscSequentialPhaseEnd(PETSC_COMM_WORLD,1);CHKERRQ(ierr);  ierr = VecView(global,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);  ierr = DMDestroy(&da);CHKERRQ(ierr);  ierr = VecDestroy(&local);CHKERRQ(ierr);  ierr = VecDestroy(&global);CHKERRQ(ierr);  ierr = PetscFinalize();  return 0;}

int main(int argc,char **argv){    Vec            u,xy;    DM             da;    PetscErrorCode ierr;    PetscInt       m = 10, n = 10, dof = 2;    PF             pf;    ierr = PetscInitialize(&argc,&argv,(char*)0,help);    CHKERRQ(ierr);    ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,m,n,PETSC_DECIDE,PETSC_DECIDE,dof,1,0,0,&da);    CHKERRQ(ierr);    ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);    CHKERRQ(ierr);    ierr = DMCreateGlobalVector(da,&u);    CHKERRQ(ierr);    ierr = DMGetCoordinates(da,&xy);    CHKERRQ(ierr);    ierr = DMDACreatePF(da,&pf);    CHKERRQ(ierr);    ierr = PFSet(pf,myfunction,0,0,0,0);    CHKERRQ(ierr);    ierr = PFSetFromOptions(pf);    CHKERRQ(ierr);    ierr = PFApplyVec(pf,xy,u);    CHKERRQ(ierr);    ierr = VecView(u,PETSC_VIEWER_DRAW_WORLD);    CHKERRQ(ierr);    /*       Free work space.  All PETSc objects should be destroyed when they       are no longer needed.    */    ierr = PFDestroy(&pf);    CHKERRQ(ierr);    ierr = DMDestroy(&da);    CHKERRQ(ierr);    ierr = PetscFinalize();    return 0;}

int main(int argc,char **argv){  TS             ts;  SNES           snes;  SNESLineSearch linesearch;  Vec            x;  AppCtx         ctx;  PetscErrorCode ierr;  DM             da;  ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;  ierr = SetFromOptions(&ctx);CHKERRQ(ierr);  ierr = TSCreate(PETSC_COMM_WORLD, &ts);CHKERRQ(ierr);  ierr = TSSetType(ts,TSCN);CHKERRQ(ierr);  ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);  ierr = TSSetIFunction(ts, NULL, FormIFunction, &ctx);CHKERRQ(ierr);  ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,N_SPECIES,1,NULL,NULL,&da);CHKERRQ(ierr);  ierr = DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);CHKERRQ(ierr);  ierr = DMDASetFieldName(da,0,"species A");CHKERRQ(ierr);  ierr = DMDASetFieldName(da,1,"species B");CHKERRQ(ierr);  ierr = DMDASetFieldName(da,2,"species C");CHKERRQ(ierr);  ierr = DMSetApplicationContext(da,&ctx);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da,&x);CHKERRQ(ierr);  ierr = FormInitialGuess(da, &ctx, x);CHKERRQ(ierr);  ierr = TSSetDM(ts, da);CHKERRQ(ierr);  ierr = TSSetDuration(ts,10000,1000.0);CHKERRQ(ierr);  ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);  ierr = TSSetInitialTimeStep(ts,0.0,1.0);CHKERRQ(ierr);  ierr = TSSetSolution(ts,x);CHKERRQ(ierr);  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);  ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);  ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr);  ierr = SNESLineSearchSetPostCheck(linesearch, ReactingFlowPostCheck, (void*)&ctx);CHKERRQ(ierr);  ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);  ierr = TSSolve(ts,x);CHKERRQ(ierr);  ierr = VecDestroy(&x);CHKERRQ(ierr);  ierr = TSDestroy(&ts);CHKERRQ(ierr);  ierr = DMDestroy(&da);CHKERRQ(ierr);  ierr = PetscFinalize();  PetscFunctionReturn(0);}

static PetscErrorCode ComputeMapping(DomainData dd,ISLocalToGlobalMapping *isg2lmap){  PetscErrorCode         ierr;  DM                     da;  AO                     ao;  DMBoundaryType         bx = DM_BOUNDARY_NONE,by = DM_BOUNDARY_NONE, bz = DM_BOUNDARY_NONE;  DMDAStencilType        stype = DMDA_STENCIL_BOX;  ISLocalToGlobalMapping temp_isg2lmap;  PetscInt               i,j,k,ig,jg,kg,lindex,gindex,localsize;  PetscInt               *global_indices;  PetscFunctionBeginUser;  /* Not an efficient mapping: this function computes a very simple lexicographic mapping     just to illustrate the creation of a MATIS object */  localsize = dd.xm_l*dd.ym_l*dd.zm_l;  ierr      = PetscMalloc1(localsize,&global_indices);CHKERRQ(ierr);  for (k=0; k<dd.zm_l; k++) {    kg=dd.startz+k;    for (j=0; j<dd.ym_l; j++) {      jg=dd.starty+j;      for (i=0; i<dd.xm_l; i++) {        ig                    =dd.startx+i;        lindex                =k*dd.xm_l*dd.ym_l+j*dd.xm_l+i;        gindex                =kg*dd.xm*dd.ym+jg*dd.xm+ig;        global_indices[lindex]=gindex;      }    }  }  if (dd.dim==3) {    ierr = DMDACreate3d(dd.gcomm,bx,by,bz,stype,dd.xm,dd.ym,dd.zm,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,NULL,&da);CHKERRQ(ierr);  } else if (dd.dim==2) {    ierr = DMDACreate2d(dd.gcomm,bx,by,stype,dd.xm,dd.ym,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);CHKERRQ(ierr);  } else {    ierr = DMDACreate1d(dd.gcomm,bx,dd.xm,1,1,NULL,&da);CHKERRQ(ierr);  }  ierr = DMDASetAOType(da,AOMEMORYSCALABLE);CHKERRQ(ierr);  ierr = DMDAGetAO(da,&ao);CHKERRQ(ierr);  ierr = AOApplicationToPetsc(ao,dd.xm_l*dd.ym_l*dd.zm_l,global_indices);CHKERRQ(ierr);  ierr = ISLocalToGlobalMappingCreate(dd.gcomm,1,localsize,global_indices,PETSC_OWN_POINTER,&temp_isg2lmap);CHKERRQ(ierr);  ierr = DMDestroy(&da);CHKERRQ(ierr);  *isg2lmap = temp_isg2lmap;  PetscFunctionReturn(0);}

PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm){  PetscInt       dim             = user->dim;  PetscBool      interpolate     = user->interpolate;  PetscReal      refinementLimit = user->refinementLimit;  const char    *partitioner     = "chaco";  PetscErrorCode ierr;  PetscFunctionBegin;  if (user->simplex) {    DM refinedMesh     = NULL;    DM distributedMesh = NULL;    ierr = DMPlexCreateBoxMesh(comm, dim, interpolate, dm);CHKERRQ(ierr);    /* Refine mesh using a volume constraint */    ierr = DMPlexSetRefinementLimit(*dm, refinementLimit);CHKERRQ(ierr);    ierr = DMRefine(*dm, comm, &refinedMesh);CHKERRQ(ierr);    if (refinedMesh) {      ierr = DMDestroy(dm);CHKERRQ(ierr);      *dm  = refinedMesh;    }    /* Distribute mesh over processes */    ierr = DMPlexDistribute(*dm, partitioner, 0, NULL, &distributedMesh);CHKERRQ(ierr);    if (distributedMesh) {      ierr = DMDestroy(dm);CHKERRQ(ierr);      *dm  = distributedMesh;    }    ierr = PetscObjectSetName((PetscObject) *dm, "Simplical Mesh");CHKERRQ(ierr);  } else {    switch (user->dim) {    case 2:      ierr = DMDACreate2d(comm, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_STENCIL_BOX, -3, -3, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, dm);CHKERRQ(ierr);      ierr = DMDASetVertexCoordinates(*dm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);CHKERRQ(ierr);      break;    default:      SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot create structured mesh of dimension %d", dim);    }    ierr = PetscObjectSetName((PetscObject) *dm, "Hexahedral Mesh");CHKERRQ(ierr);  }  ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);  PetscFunctionReturn(0);}

int main(int argc,char **argv){  TS                     ts;  Vec                    x,c;  PetscErrorCode         ierr;  DM                     da;  PetscInitialize(&argc,&argv,(char *)0,help);  ierr = TSCreate(PETSC_COMM_WORLD, &ts);CHKERRQ(ierr);  ierr = TSSetType(ts,TSCN);CHKERRQ(ierr);  ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);  ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,&da);CHKERRQ(ierr);  ierr = DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);CHKERRQ(ierr);  ierr = DMDASetFieldName(da,0,"Heat");CHKERRQ(ierr);  ierr = DMCreateGlobalVector(da,&x);CHKERRQ(ierr);  ierr = FormInitialGuess(da,PETSC_NULL,x);CHKERRQ(ierr);  ierr = DMDATSSetIFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,PetscReal,void*,void*,void*,void*))FormIFunctionLocal,PETSC_NULL);CHKERRQ(ierr);  /* set up the coefficient */  ierr = DMGetNamedGlobalVector(da,"coefficient",&c);CHKERRQ(ierr);  ierr = FormDiffusionCoefficient(da,PETSC_NULL,c);CHKERRQ(ierr);  ierr = DMRestoreNamedGlobalVector(da,"coefficient",&c);CHKERRQ(ierr);  ierr = DMCoarsenHookAdd(da,PETSC_NULL,CoefficientRestrictHook,ts);CHKERRQ(ierr);  ierr = DMSubDomainHookAdd(da,PETSC_NULL,CoefficientSubDomainRestrictHook,ts);CHKERRQ(ierr);  ierr = TSSetDM(ts, da);CHKERRQ(ierr);  ierr = TSSetDuration(ts,10000,1000.0);CHKERRQ(ierr);  ierr = TSSetInitialTimeStep(ts,0.0,0.05);CHKERRQ(ierr);  ierr = TSSetSolution(ts,x);CHKERRQ(ierr);  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);  ierr = TSSolve(ts,x);CHKERRQ(ierr);  ierr = VecDestroy(&x);CHKERRQ(ierr);  ierr = TSDestroy(&ts);CHKERRQ(ierr);  ierr = DMDestroy(&da);CHKERRQ(ierr);  ierr = PetscFinalize();  PetscFunctionReturn(0);}

int main(int argc,char **argv){  KSP            ksp;  DM             da;  Vec            x, xNew;  UserContext    user;  PetscErrorCode ierr;  PetscInitialize(&argc,&argv,(char*)0,help);  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);  ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,3,3,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);CHKERRQ(ierr);  ierr = DMSetApplicationContext(da, &user);CHKERRQ(ierr);  ierr = KSPSetDM(ksp, da);CHKERRQ(ierr);  ierr     = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for PCICE", "DM");  user.phi = 0.5;  ierr     = PetscOptionsScalar("-phi", "The time weighting parameter", "ex31.c", user.phi, &user.phi, NULL);CHKERRQ(ierr);  user.dt  = 0.1;  ierr     = PetscOptionsScalar("-dt", "The time step", "ex31.c", user.dt, &user.dt, NULL);CHKERRQ(ierr);  ierr     = PetscOptionsEnd();  ierr = CreateStructures(da, &user);CHKERRQ(ierr);  ierr = ComputePredictor(da, &user);CHKERRQ(ierr);  ierr = KSPSetComputeRHS(ksp,ComputeRHS,&user);CHKERRQ(ierr);  ierr = KSPSetComputeOperators(ksp,ComputeMatrix,&user);CHKERRQ(ierr);  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);  ierr = KSPSolve(ksp, NULL, NULL);CHKERRQ(ierr);  ierr = KSPGetSolution(ksp, &x);CHKERRQ(ierr);  ierr = VecDuplicate(x, &xNew);CHKERRQ(ierr);  ierr = ComputeCorrector(da, x, xNew);CHKERRQ(ierr);  ierr = VecDestroy(&xNew);CHKERRQ(ierr);  ierr = DestroyStructures(da, &user);CHKERRQ(ierr);  ierr = DMDestroy(&da);CHKERRQ(ierr);  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);  ierr = PetscFinalize();  return 0;}

PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm){  PetscInt       dim = user->dim;  PetscErrorCode ierr;  PetscFunctionBegin;  switch(dim) {  case 2:    ierr = DMDACreate2d(comm, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_STENCIL_STAR, -4, -4, PETSC_DECIDE, PETSC_DECIDE, 2, 1, PETSC_NULL, PETSC_NULL, dm);CHKERRQ(ierr);    break;  case 3:    ierr = DMDACreate3d(comm, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_STENCIL_STAR, -4, -4, -4, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 2, 1, PETSC_NULL, PETSC_NULL, PETSC_NULL, dm);CHKERRQ(ierr);    break;  default:    SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Could not create mesh for dimension %d", dim);  }  ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);  ierr = DMDASetFieldName(*dm, 0, "velocity");CHKERRQ(ierr);  ierr = DMDASetFieldName(*dm, 1, "pressure");CHKERRQ(ierr);  user->dm = *dm;  PetscFunctionReturn(0);}