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

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

PetscErrorCode DADefineXLinearField2D(DM da,Vec field){  PetscErrorCode ierr;  PetscInt       i,j;  PetscInt       sx,nx,sy,ny;  Vec            Gcoords;  DMDACoor2d     **XX;  PetscScalar    **FF;  DM             cda;  PetscFunctionBeginUser;  ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr);  ierr = DMGetCoordinates(da,&Gcoords);CHKERRQ(ierr);  ierr = DMDAVecGetArray(cda,Gcoords,&XX);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,field,&FF);CHKERRQ(ierr);  ierr = DMDAGetCorners(da,&sx,&sy,0,&nx,&ny,0);CHKERRQ(ierr);  for (i=sx; i<sx+nx; i++) {    for (j=sy; j<sy+ny; j++ ) {      FF[j][i] = 10.0 + 3.0 * XX[j][i].x + 5.5 * XX[j][i].y + 8.003 * XX[j][i].x * XX[j][i].y;    }  }  ierr = DMDAVecRestoreArray(da,field,&FF);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(cda,Gcoords,&XX);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode ReactingFlowPostCheck(SNESLineSearch linesearch, Vec X, Vec Y, Vec W, PetscBool *changed_y, PetscBool *changed_w, void *vctx){  PetscInt       i,j,l,Mx,My,xs,ys,xm,ym;  PetscErrorCode ierr;  Field          **x;  SNES           snes;  DM             da;  PetscScalar    min;  PetscFunctionBeginUser;   *changed_w = PETSC_FALSE;  ierr = VecMin(X,NULL,&min);CHKERRQ(ierr);  if (min >= 0.) PetscFunctionReturn(0);  *changed_w = PETSC_TRUE;  ierr = SNESLineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr);  ierr = SNESGetDM(snes,&da);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);  ierr = DMDAVecGetArray(da,W,&x);CHKERRQ(ierr);  ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);  for (j=ys; j<ys+ym; j++) {    for (i=xs; i<xs+xm; i++) {      for (l = 0; l < N_SPECIES; l++) {        if (x[j][i].sp[l] < 0.) x[j][i].sp[l] = 0.;      }    }  }  ierr = DMDAVecRestoreArray(da,W,&x);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode FillMatrix(DM da,Mat A){  PetscErrorCode ierr;  PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,idx;  PetscScalar    v[7];  MatStencil     row,col[7];  PetscFunctionBeginUser;  ierr = DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);  ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr);  for (k=zs;k<zs+zm;k++) {    for (j=ys;j<ys+ym;j++) {      for (i=xs;i<xs+xm;i++) {        row.i=i; row.j=j; row.k=k;        col[0].i=row.i; col[0].j=row.j; col[0].k=row.k;        v[0]=6.0;        idx=1;        if (k>0) { v[idx]=-1.0; col[idx].i=i; col[idx].j=j; col[idx].k=k-1; idx++; }        if (j>0) { v[idx]=-1.0; col[idx].i=i; col[idx].j=j-1; col[idx].k=k; idx++; }        if (i>0) { v[idx]=-1.0; col[idx].i=i-1; col[idx].j=j; col[idx].k=k; idx++; }        if (i<mx-1) { v[idx]=-1.0; col[idx].i=i+1; col[idx].j=j; col[idx].k=k; idx++; }        if (j<my-1) { v[idx]=-1.0; col[idx].i=i; col[idx].j=j+1; col[idx].k=k; idx++; }        if (k<mz-1) { v[idx]=-1.0; col[idx].i=i; col[idx].j=j; col[idx].k=k+1; idx++; }        ierr = MatSetValuesStencil(A,1,&row,idx,col,v,INSERT_VALUES);CHKERRQ(ierr);      }    }  }  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode ComputeB(AppCtx *user){  PetscErrorCode ierr;  PetscInt       i,j;  PetscInt       nx,ny,xs,xm,ys,ym;  PetscReal      two=2.0, pi=4.0*atan(1.0);  PetscReal      hx,hy,ehxhy;  PetscReal      temp;  PetscReal      ecc=user->ecc;  PetscReal      **b;  PetscFunctionBeginUser;  nx    = user->nx;  ny    = user->ny;  hx    = two*pi/(nx+1.0);  hy    = two*user->b/(ny+1.0);  ehxhy = ecc*hx*hy;  /* Get pointer to local vector data */  ierr = DMDAVecGetArray(user->da,user->B, &b);CHKERRQ(ierr);  ierr = DMDAGetCorners(user->da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);  /* Compute the linear term in the objective function */  for (i=xs; i<xs+xm; i++) {    temp=PetscSinReal((i+1)*hx);    for (j=ys; j<ys+ym; j++) b[j][i] = -ehxhy*temp;  }  /* Restore vectors */  ierr = DMDAVecRestoreArray(user->da,user->B,&b);CHKERRQ(ierr);  ierr = PetscLogFlops(5*xm*ym+3*xm);CHKERRQ(ierr);  PetscFunctionReturn(0);}

static PetscErrorCode ComputeMatrix(KSP ksp,Mat J,Mat jac,MatStructure *str,void *ctx){  AppCtx         *user = (AppCtx*)ctx;  PetscErrorCode ierr;  PetscInt       i,mx,xm,xs;  PetscScalar    v[3],h,xlow,xhigh;  MatStencil     row,col[3];  DM             da;  PetscFunctionBeginUser;  ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);  ierr = DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);  ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr);  h    = 1.0/(mx-1);  for (i=xs; i<xs+xm; i++){    row.i = i;    if (i==0 || i==mx-1){      v[0] = 2.0;      ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);    } else {       xlow  = h*(PetscReal)i - .5*h;       xhigh = xlow + h;       v[0] = (-1.0 - user->e*PetscSinScalar(2.0*PETSC_PI*user->k*xlow))/h;col[0].i = i-1;       v[1] = (2.0 + user->e*PetscSinScalar(2.0*PETSC_PI*user->k*xlow) + user->e*PetscSinScalar(2.0*PETSC_PI*user->k*xhigh))/h;col[1].i = row.i;       v[2] = (-1.0 - user->e*PetscSinScalar(2.0*PETSC_PI*user->k*xhigh))/h;col[2].i = i+1;      ierr = MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);CHKERRQ(ierr);    }  }  ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx){  PetscErrorCode ierr;  PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;  DM             dm;  PetscScalar    Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy;  PetscScalar    ***barray;  PetscFunctionBeginUser;  ierr    = KSPGetDM(ksp,&dm);CHKERRQ(ierr);  ierr    = DMDAGetInfo(dm,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);  Hx      = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);  HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;  ierr    = DMDAGetCorners(dm,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr);  ierr    = DMDAVecGetArray(dm,b,&barray);CHKERRQ(ierr);  for (k=zs; k<zs+zm; k++) {    for (j=ys; j<ys+ym; j++) {      for (i=xs; i<xs+xm; i++) {        if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) {          barray[k][j][i] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);        } else {          barray[k][j][i] = Hx*Hy*Hz;        }      }    }  }  ierr = DMDAVecRestoreArray(dm,b,&barray);CHKERRQ(ierr);  PetscFunctionReturn(0);}

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

PetscErrorCode FormInitialGuess(DM da,AppCtx *ctx,Vec X){  PetscInt       i,j,l,Mx,My,xs,ys,xm,ym;  PetscErrorCode ierr;  Field          **x;  PetscFunctionBeginUser;  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);  ierr = DMDAVecGetArray(da,X,&x);CHKERRQ(ierr);  ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);  for (j=ys; j<ys+ym; j++) {    for (i=xs; i<xs+xm; i++) {      for (l = 0; l < N_SPECIES; l++) {        if (i == 0) {          if (l == 0)      x[j][i].sp[l] = (ctx->x_inflow.sp[l]*((PetscScalar)j) / (My - 1));          else if (l == 1) x[j][i].sp[l] = (ctx->x_inflow.sp[l]*(1. - ((PetscScalar)j) / (My - 1)));          else             x[j][i].sp[l] = ctx->x_0.sp[l];        }      }    }  }  ierr = DMDAVecRestoreArray(da,X,&x);CHKERRQ(ierr);  PetscFunctionReturn(0);}

/*    FormInitialGuess - Forms initial approximation.   Input Parameters:   user - user-defined application context   X - vector   Output Parameter:   X - vector */PetscErrorCode FormInitialGuess(AppCtx *user,Vec X){  PetscInt       i,j,row,mx,my,xs,ys,xm,ym,gxm,gym,gxs,gys;  PetscErrorCode ierr;  PetscReal      one = 1.0,lambda,temp1,temp,hx,hy,hxdhy,hydhx,sc;  PetscScalar    *x;  Vec            localX = user->localX;  mx = user->mx;            my = user->my;            lambda = user->param;  hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);  sc = hx*hy*lambda;        hxdhy = hx/hy;            hydhx = hy/hx;  temp1 = lambda/(lambda + one);  /*     Get a pointer to vector data.       - For default PETSc vectors,VecGetArray() returns a pointer to         the data array.  Otherwise, the routine is implementation dependent.       - You MUST call VecRestoreArray() when you no longer need access to         the array.  */  ierr = VecGetArray(localX,&x);CHKERRQ(ierr);  /*     Get local grid boundaries (for 2-dimensional DMDA):       xs, ys   - starting grid indices (no ghost points)       xm, ym   - widths of local grid (no ghost points)       gxs, gys - starting grid indices (including ghost points)       gxm, gym - widths of local grid (including ghost points)  */  ierr = DMDAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);CHKERRQ(ierr);  ierr = DMDAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);CHKERRQ(ierr);  /*     Compute initial guess over the locally owned part of the grid  */  for (j=ys; j<ys+ym; j++) {    temp = (PetscReal)(PetscMin(j,my-j-1))*hy;    for (i=xs; i<xs+xm; i++) {      row = i - gxs + (j - gys)*gxm;       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {        x[row] = 0.0;         continue;      }      x[row] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));     }  }  /*     Restore vector  */  ierr = VecRestoreArray(localX,&x);CHKERRQ(ierr);  /*     Insert values into global vector  */  ierr = DMLocalToGlobalBegin(user->da,localX,INSERT_VALUES,X);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(user->da,localX,INSERT_VALUES,X);CHKERRQ(ierr);  return 0;} 

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

PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat *J,Mat *Jpre,MatStructure *str,void *ctx){  PetscErrorCode ierr;  PetscInt       i,j,Mx,My,xs,ys,xm,ym,nc;  AppCtx         *user = (AppCtx*)ctx;  DM             da = (DM)user->da;  MatStencil     col[5],row;  PetscScalar    vals[5],hx,hy,sx,sy;  PetscFunctionBeginUser;  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);  ierr = DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);CHKERRQ(ierr);  hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);  hy = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy);  for (j=ys; j<ys+ym; j++){    for (i=xs; i<xs+xm; i++){      nc    = 0;      row.j = j; row.i = i;      if (user->boundary == 0 && (i == 0 || i == Mx-1 || j == 0 || j == My-1)) {        col[nc].j = j; col[nc].i = i; vals[nc++] = 1.0;      } else if (user->boundary > 0 && i == 0) {  /* Left Neumann */        col[nc].j = j; col[nc].i = i;   vals[nc++] = 1.0;        col[nc].j = j; col[nc].i = i+1; vals[nc++] = -1.0;      } else if (user->boundary > 0 && i == Mx-1){/* Right Neumann */        col[nc].j = j; col[nc].i = i;   vals[nc++] = 1.0;        col[nc].j = j; col[nc].i = i-1; vals[nc++] = -1.0;      } else if (user->boundary > 0 && j == 0) {  /* Bottom Neumann */        col[nc].j = j;   col[nc].i = i; vals[nc++] = 1.0;        col[nc].j = j+1; col[nc].i = i; vals[nc++] = -1.0;      } else if (user->boundary > 0 && j == My-1){/* Top Neumann */        col[nc].j = j;   col[nc].i = i;  vals[nc++] = 1.0;        col[nc].j = j-1; col[nc].i = i;  vals[nc++] = -1.0;      } else {   /* Interior */        col[nc].j = j-1; col[nc].i = i;   vals[nc++] = -sy;        col[nc].j = j;   col[nc].i = i-1; vals[nc++] = -sx;        col[nc].j = j;   col[nc].i = i;   vals[nc++] = 2.0*(sx + sy) + a;        col[nc].j = j;   col[nc].i = i+1; vals[nc++] = -sx;        col[nc].j = j+1; col[nc].i = i;   vals[nc++] = -sy;      }      ierr = MatSetValuesStencil(*Jpre,1,&row,nc,col,vals,INSERT_VALUES);CHKERRQ(ierr);    }  }  ierr = MatAssemblyBegin(*Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  ierr = MatAssemblyEnd(*Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  if (*J != *Jpre) {    ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);    ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  }  if (user->viewJacobian){    ierr = PetscPrintf(((PetscObject)*Jpre)->comm,"Jpre:/n");CHKERRQ(ierr);    ierr = MatView(*Jpre,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);  }  PetscFunctionReturn(0);}

/*  Compute the right thickness H = H(x).  See lecture notes for exact solution. */PetscErrorCode FillThicknessAndExactSoln(DM da, AppCtx *user){  PetscErrorCode ierr;  PetscInt       i,Mx,xs,xm;  PetscReal      hx, n, r, Cs, xx, flux, qg, p, B, dudx, ustag;  PetscScalar    *H, *uex, *visc;  PetscFunctionBegin;  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);  ierr = DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);  /* constants, independent of x */  hx = user->L / (PetscReal)(Mx-1);  n  = user->n;  r  = user->rho / user->rhow;  Cs = user->A * PetscPowScalar( ( 0.25 * user->rho * user->g * (1.0 - r) ), n);  qg = user->ug * user->Hg;  p  = 1.0 + 1.0 / user->n;  B  = PetscPowScalar(user->A,-1.0/user->n);  /* Compute regular grid exact soln and staggered-grid thickness over the     locally-owned part of the grid */  ierr = DMDAVecGetArray(da,user->uexact,&uex);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,user->H,&H);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da,user->viscosity,&visc);CHKERRQ(ierr);  for (i=xs; i<xs+xm; i++) {    /* get exact velocity and strain rate on regular grid */    xx = hx * (PetscReal)i;  /* = x_i = distance from grounding line */    flux = user->accum * xx + qg; /* flux at x_i */    GetUEx(user->ug, qg, user->accum, n, Cs, flux, &(uex[i]), &dudx);    /* exact viscosity on regular grid */    visc[i] = GetViscosityFromStrainRate(dudx, p, B);    /* exact thickness on staggered grid */    flux += user->accum * hx * 0.5; /* flux at x_{i+1/2} */    GetUEx(user->ug, qg, user->accum, n, Cs, flux, &ustag, &dudx);    H[i] = flux / ustag;  }  ierr = DMDAVecRestoreArray(da,user->uexact,&uex);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,user->H,&H);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da,user->viscosity,&visc);CHKERRQ(ierr);  /* separately compute and store calving-front values */  flux = user->accum * user->L + qg;  GetUEx(user->ug, qg, user->accum, n, Cs, flux, &(user->uexactcalv), &dudx);  user->Hcalv = flux / user->uexactcalv;  /* strain rate at calving front */  /* MATLAB: gamma = ( 0.25 * p.A^(1/n) * (1 - r) * p.rho * p.g * H(end) )^n; */  user->gamma = 0.25 * PetscPowScalar(user->A,1.0/user->n) * (1.0 - r)                  * user->rho * user->g * user->Hcalv;  user->gamma = PetscPowScalar(user->gamma,user->n);  PetscFunctionReturn(0);}

/*   RHSMatrixLaplacian - User-provided routine to compute the right-hand-side   matrix for the Laplacian operator   Input Parameters:   ts - the TS context   t - current time  (ignored)   X - current solution (ignored)   dummy - optional user-defined context, as set by TSetRHSJacobian()   Output Parameters:   AA - Jacobian matrix   BB - optionally different matrix from which the preconditioner is built   str - flag indicating matrix structure*/PetscErrorCode RHSMatrixLaplaciangllDM(TS ts,PetscReal t,Vec X,Mat A,Mat BB,void *ctx){  PetscReal      **temp;  PetscReal      vv;  AppCtx         *appctx = (AppCtx*)ctx;     /* user-defined application context */  PetscErrorCode ierr;  PetscInt       i,xs,xn,l,j;  PetscInt       *rowsDM;  PetscBool      flg = PETSC_FALSE;  ierr = PetscOptionsGetBool(NULL,NULL,"-gll_mf",&flg,NULL);CHKERRQ(ierr);  if (!flg) {    /*     Creates the element stiffness matrix for the given gll     */    ierr = PetscGLLElementLaplacianCreate(&appctx->SEMop.gll,&temp);CHKERRQ(ierr);    /* workarround for clang analyzer warning: Division by zero */    if (appctx->param.N <= 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Spectral element order should be > 1");    /* scale by the size of the element */    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 = MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);    ierr = DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);CHKERRQ(ierr);    xs   = xs/(appctx->param.N-1);    xn   = xn/(appctx->param.N-1);    ierr = PetscMalloc1(appctx->param.N,&rowsDM);CHKERRQ(ierr);    /*     loop over local elements     */    for (j=xs; j<xs+xn; j++) {      for (l=0; l<appctx->param.N; l++) {        rowsDM[l] = 1+(j-xs)*(appctx->param.N-1)+l;      }      ierr = MatSetValuesLocal(A,appctx->param.N,rowsDM,appctx->param.N,rowsDM,&temp[0][0],ADD_VALUES);CHKERRQ(ierr);    }    ierr = PetscFree(rowsDM);CHKERRQ(ierr);    ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);    ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);    ierr = VecReciprocal(appctx->SEMop.mass);CHKERRQ(ierr);    ierr = MatDiagonalScale(A,appctx->SEMop.mass,0);CHKERRQ(ierr);    ierr = VecReciprocal(appctx->SEMop.mass);CHKERRQ(ierr);    ierr = PetscGLLElementLaplacianDestroy(&appctx->SEMop.gll,&temp);CHKERRQ(ierr);  } else {    ierr = MatSetType(A,MATSHELL);CHKERRQ(ierr);    ierr = MatSetUp(A);CHKERRQ(ierr);    ierr = MatShellSetContext(A,appctx);CHKERRQ(ierr);    ierr = MatShellSetOperation(A,MATOP_MULT,(void (*)(void))MatMult_Laplacian);CHKERRQ(ierr);  }  return 0;}

/*   FormInitialGuess - Forms initial approximation.   Input Parameters:   user - user-defined application context   X - vector   Output Parameter:   X - vector */PetscErrorCode FormInitialGuess(AppCtx *user,Vec X){  PetscInt       i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm;  PetscErrorCode ierr;  PetscReal      lambda,temp1,hx,hy,hz,tempk,tempj;  PetscScalar    ***x;  PetscFunctionBeginUser;  ierr = DMDAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,&Mz,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);  hz     = 1.0/(PetscReal)(Mz-1);  temp1  = lambda/(lambda + 1.0);  /*     Get a pointer to vector data.       - For default PETSc vectors, VecGetArray() returns a pointer to         the data array.  Otherwise, the routine is implementation dependent.       - You MUST call VecRestoreArray() when you no longer need access to         the array.  */  ierr = DMDAVecGetArray(user->da,X,&x);CHKERRQ(ierr);  /*     Get local grid boundaries (for 3-dimensional DMDA):       xs, ys, zs   - starting grid indices (no ghost points)       xm, ym, zm   - widths of local grid (no ghost points)  */  ierr = DMDAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr);  /*     Compute initial guess over the locally owned part of the grid  */  for (k=zs; k<zs+zm; k++) {    tempk = (PetscReal)(PetscMin(k,Mz-k-1))*hz;    for (j=ys; j<ys+ym; j++) {      tempj = PetscMin((PetscReal)(PetscMin(j,My-j-1))*hy,tempk);      for (i=xs; i<xs+xm; i++) {        if (i == 0 || j == 0 || k == 0 || i == Mx-1 || j == My-1 || k == Mz-1) {          /* boundary conditions are all zero Dirichlet */          x[k][j][i] = 0.0;        } else {          x[k][j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,tempj));        }      }    }  }  /*     Restore vector  */  ierr = DMDAVecRestoreArray(user->da,X,&x);CHKERRQ(ierr);  PetscFunctionReturn(0);}

	PetscErrorCode TcSolver::generateBC2(){	PetscErrorCode	ierr;	PetscInt	mstart, nstart, m, n, i, j;	PetscReal	**qx, **qy;	PetscReal	**bc22;	Vec		bc2Global;	ierr = VecSet(bc2, 0.0); CHKERRQ(ierr);	ierr = DMCompositeGetAccess(lambdaPack, bc2, &bc2Global, NULL); CHKERRQ(ierr);		ierr = DMDAVecGetArray(uda, qxLocal, &qx); CHKERRQ(ierr);	ierr = DMDAVecGetArray(vda, qyLocal, &qy); CHKERRQ(ierr);	ierr = DMDAVecGetArray(pda, bc2Global, &bc22); CHKERRQ(ierr);	ierr = DMDAGetCorners(pda, &mstart, &nstart, NULL, &m, &n, NULL); CHKERRQ(ierr);	//x-faces		for(j=nstart; j<nstart+n; j++)	{			//-X		if(mstart == 0) //if -X is current process		{			bc22[j][0] -= qx[j][-1];		}		//+X		if(mstart+m-1 == fluid.nx-1)		{			bc22[j][fluid.nx-1] += qx[j][fluid.nx-1];		}	}		//y-faces	for(i=mstart; i<mstart+m; i++)	{		//-Y		if(nstart ==0)		{			bc22[0][i] -= qy[-1][i];		}		//+Y		if(nstart+n-1 == fluid.ny-1)		{			bc22[fluid.ny-1][i] += qy[fluid.ny-1][i];		}	}		ierr = DMDAVecRestoreArray(pda, bc2Global, &bc22); CHKERRQ(ierr);	ierr = DMDAVecRestoreArray(uda, qxLocal, &qx); CHKERRQ(ierr);	ierr = DMDAVecRestoreArray(vda, qyLocal, &qy); CHKERRQ(ierr);	ierr = DMCompositeRestoreAccess(lambdaPack, bc2, &bc2Global, NULL); CHKERRQ(ierr);	return 0;}

/*      Evaluates FU = Gradiant(L(w,u,lambda))*/PetscErrorCode FormFunction(SNES snes,Vec U,Vec FU,void *dummy){  UserCtx        *user = (UserCtx*)dummy;  PetscErrorCode ierr;  PetscInt       xs,xm,i,N;  PetscScalar    *u,*lambda,*w,*fu,*fw,*flambda,d,h;  Vec            vw,vu,vlambda,vfw,vfu,vflambda;  PetscFunctionBeginUser;  ierr = DMCompositeGetLocalVectors(user->packer,&vw,&vu,&vlambda);CHKERRQ(ierr);  ierr = DMCompositeGetLocalVectors(user->packer,&vfw,&vfu,&vflambda);CHKERRQ(ierr);  ierr = DMCompositeScatter(user->packer,U,vw,vu,vlambda);CHKERRQ(ierr);  ierr = DMDAGetCorners(user->da1,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);  ierr = DMDAGetInfo(user->da1,0,&N,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);  ierr = VecGetArray(vw,&w);CHKERRQ(ierr);  ierr = VecGetArray(vfw,&fw);CHKERRQ(ierr);  ierr = DMDAVecGetArray(user->da1,vu,&u);CHKERRQ(ierr);  ierr = DMDAVecGetArray(user->da1,vfu,&fu);CHKERRQ(ierr);  ierr = DMDAVecGetArray(user->da1,vlambda,&lambda);CHKERRQ(ierr);  ierr = DMDAVecGetArray(user->da1,vflambda,&flambda);CHKERRQ(ierr);  d    = (N-1.0);  h    = 1.0/d;  /* derivative of L() w.r.t. w */  if (xs == 0) { /* only first processor computes this */    fw[0] = -2.*d*lambda[0];  }  /* derivative of L() w.r.t. u */  for (i=xs; i<xs+xm; i++) {    if      (i == 0)   flambda[0]   =    h*u[0]   + 2.*d*lambda[0]   - d*lambda[1];    else if (i == 1)   flambda[1]   = 2.*h*u[1]   + 2.*d*lambda[1]   - d*lambda[2];    else if (i == N-1) flambda[N-1] =    h*u[N-1] + 2.*d*lambda[N-1] - d*lambda[N-2];    else if (i == N-2) flambda[N-2] = 2.*h*u[N-2] + 2.*d*lambda[N-2] - d*lambda[N-3];    else               flambda[i]   = 2.*h*u[i]   - d*(lambda[i+1] - 2.0*lambda[i] + lambda[i-1]);  }  /* derivative of L() w.r.t. lambda */  for (i=xs; i<xs+xm; i++) {    if      (i == 0)   fu[0]   = 2.0*d*(u[0] - w[0]);    else if (i == N-1) fu[N-1] = 2.0*d*u[N-1];    else               fu[i]   = -(d*(u[i+1] - 2.0*u[i] + u[i-1]) - 2.0*h);  }  ierr = VecRestoreArray(vw,&w);CHKERRQ(ierr);  ierr = VecRestoreArray(vfw,&fw);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(user->da1,vu,&u);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(user->da1,vfu,&fu);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(user->da1,vlambda,&lambda);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(user->da1,vflambda,&flambda);CHKERRQ(ierr);  ierr = DMCompositeGather(user->packer,FU,INSERT_VALUES,vfw,vfu,vflambda);CHKERRQ(ierr);  ierr = DMCompositeRestoreLocalVectors(user->packer,&vw,&vu,&vlambda);CHKERRQ(ierr);  ierr = DMCompositeRestoreLocalVectors(user->packer,&vfw,&vfu,&vflambda);CHKERRQ(ierr);  PetscFunctionReturn(0);}

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 efs_set_sol(efs *slv, double sol[]){	PetscErrorCode ierr;	ierr = DMCreateGlobalVector(slv->dm, &slv->sol);	slv->state.sol = sol;	if (slv->grid.nd == 2) {		PetscScalar **svec;		double **sol_state;		PetscInt i, j;		PetscInt xs, ys, xm, ym;		ierr = DMDAGetCorners(slv->dm, &xs, &ys, 0, &xm, &ym, 0);CHKERRQ(ierr);		ierr = ef_dmap_get(slv->dmap, slv->state.sol, &sol_state);CHKERRQ(ierr);		ierr = DMDAVecGetArray(slv->dm, slv->sol, &svec);CHKERRQ(ierr);		for (j = ys; j < ys + ym; j++) {			for (i = xs; i < xs + xm; i++) {				svec[j][i] = sol_state[j][i];			}		}		ierr = ef_dmap_restore(slv->dmap, &sol_state);CHKERRQ(ierr);		ierr = DMDAVecRestoreArray(slv->dm, slv->sol, &svec);CHKERRQ(ierr);	} else if (slv->grid.nd == 3) {		PetscScalar ***svec;		double ***sol_state;		PetscInt i, j, k;		PetscInt xs, ys, zs, xm, ym, zm;		ierr = DMDAGetCorners(slv->dm, &xs, &ys, &zs, &xm, &ym, &zm);CHKERRQ(ierr);		ierr = ef_dmap_get(slv->dmap, slv->state.sol, &sol_state);CHKERRQ(ierr);		ierr = DMDAVecGetArray(slv->dm, slv->sol, &svec);CHKERRQ(ierr);		for (k = zs; k < zs + zm; k++) {			for (j = ys; j < ys + ym; j++) {				for (i = xs; i < xs + xm; i++) {					svec[k][j][i] = sol_state[k][j][i];				}			}		}		ierr = ef_dmap_restore(slv->dmap, &sol_state);CHKERRQ(ierr);		ierr = DMDAVecRestoreArray(slv->dm, slv->sol, &svec);CHKERRQ(ierr);	}	return 0;}

/*  This solves the  equation Δu = qρ with  Laplacian defined by 7-point  stencil and  periodic boundary  conditions.  Thus with  q = 1  it is  exactly the pseudo-inverse of  the Laplacian operator implemented by  7-point  stencil.   Given  the  periodicity the  pseudo-inverse  may  differ by at most a constant.  To get the potential in kcal/mol as used in the rest of the code you  need to supply q = -4π/ε
C++ DMDAVecRestoreArray函数代码示例
C++ DMDACreate2d函数代码示例