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

PetscErrorCode DMView_DA_VTK(DM da, PetscViewer viewer){  PetscInt       dim, dof, M = 0, N = 0, P = 0;  PetscErrorCode ierr;  PetscFunctionBegin;  ierr = DMDAGetInfo(da, &dim, &M, &N, &P, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);CHKERRQ(ierr);  if (!da->coordinates) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP, "VTK output requires DMDA coordinates.");  /* Write Header */  ierr = PetscViewerASCIIPrintf(viewer,"# vtk DataFile Version 2.0/n");CHKERRQ(ierr);  ierr = PetscViewerASCIIPrintf(viewer,"Structured Mesh Example/n");CHKERRQ(ierr);  ierr = PetscViewerASCIIPrintf(viewer,"ASCII/n");CHKERRQ(ierr);  ierr = PetscViewerASCIIPrintf(viewer,"DATASET STRUCTURED_GRID/n");CHKERRQ(ierr);  ierr = PetscViewerASCIIPrintf(viewer,"DIMENSIONS %d %d %d/n", M, N, P);CHKERRQ(ierr);  ierr = PetscViewerASCIIPrintf(viewer,"POINTS %d double/n", M*N*P);CHKERRQ(ierr);  if (da->coordinates) {    DM  dac;    Vec natural;    ierr = DMGetCoordinateDM(da, &dac);CHKERRQ(ierr);    ierr = DMDACreateNaturalVector(dac, &natural);CHKERRQ(ierr);    ierr = PetscObjectSetOptionsPrefix((PetscObject) natural, "coor_");CHKERRQ(ierr);    ierr = DMDAGlobalToNaturalBegin(dac, da->coordinates, INSERT_VALUES, natural);CHKERRQ(ierr);    ierr = DMDAGlobalToNaturalEnd(dac, da->coordinates, INSERT_VALUES, natural);CHKERRQ(ierr);    ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_VTK_COORDS);CHKERRQ(ierr);    ierr = VecView(natural, viewer);CHKERRQ(ierr);    ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);    ierr = VecDestroy(&natural);CHKERRQ(ierr);  }  PetscFunctionReturn(0);}

PetscErrorCode FormPsiAndInitialGuess(DM da,Vec U0,PetscBool feasible){  ObsCtx         *user;  PetscErrorCode ierr;  PetscInt       i,j,Mx,My,xs,ys,xm,ym;  DM             coordDA;  Vec            coordinates;  DMDACoor2d     **coords;  PetscReal      **psi, **u0, **uexact,                 x, y, r,                 afree = 0.69797, A = 0.68026, B = 0.47152, pi = 3.1415926;  PetscFunctionBeginUser;  ierr = DMGetApplicationContext(da,&user);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 = DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);CHKERRQ(ierr);  ierr = DMGetCoordinateDM(da, &coordDA);CHKERRQ(ierr);  ierr = DMGetCoordinates(da, &coordinates);CHKERRQ(ierr);  ierr = DMDAVecGetArray(coordDA, coordinates, &coords);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da, user->psi, &psi);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da, U0, &u0);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da, user->uexact, &uexact);CHKERRQ(ierr);  for (j=ys; j<ys+ym; j++) {    for (i=xs; i<xs+xm; i++) {      x = coords[j][i].x;      y = coords[j][i].y;      r = sqrt(x * x + y * y);      if (r <= 1.0) {        psi[j][i] = sqrt(1.0 - r * r);      } else {        psi[j][i] = -1.0;      }      if (r <= afree) {        uexact[j][i] = psi[j][i];  /* on the obstacle */      } else {        uexact[j][i] = - A * log(r) + B;   /* solves the laplace eqn */      }      if (feasible) {        if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {          u0[j][i] = uexact[j][i];        } else {          /* initial guess is admissible: it is above the obstacle */          u0[j][i] = uexact[j][i] + cos(pi * x / 4) * cos(pi * y / 4);        }      } else {        u0[j][i] = 0.;      }    }  }  ierr = DMDAVecRestoreArray(da, user->psi, &psi);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da, U0, &u0);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da, user->uexact, &uexact);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(coordDA, coordinates, &coords);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode FormCoordinates(DM da,AppCtx *user){  PetscErrorCode ierr;  Vec            coords;  DM             cda;  PetscInt       mx,my,mz;  PetscInt       i,j,k,xs,ys,zs,xm,ym,zm;  CoordField     ***x;  PetscFunctionBegin;  ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr);  ierr = DMCreateGlobalVector(cda,&coords);CHKERRQ(ierr);  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);  ierr = DMDAVecGetArray(da,coords,&x);CHKERRQ(ierr);  for (k=zs; k<zs+zm; k++) {    for (j=ys; j<ys+ym; j++) {      for (i=xs; i<xs+xm; i++) {        PetscReal cx = ((PetscReal)i) / (((PetscReal)(mx-1)));        PetscReal cy = ((PetscReal)j) / (((PetscReal)(my-1)));        PetscReal cz = ((PetscReal)k) / (((PetscReal)(mz-1)));        PetscReal rad = user->rad + cy*user->height;        PetscReal ang = (cx - 0.5)*user->arc;        x[k][j][i][0] = rad*PetscSinReal(ang);        x[k][j][i][1] = rad*PetscCosReal(ang) - (user->rad + 0.5*user->height)*PetscCosReal(-0.5*user->arc);        x[k][j][i][2] = user->width*(cz - 0.5);      }    }  }  ierr = DMDAVecRestoreArray(da,coords,&x);CHKERRQ(ierr);  ierr = DMSetCoordinates(da,coords);CHKERRQ(ierr);  ierr = VecDestroy(&coords);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode DMDAComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal detJ[]){  DM               cdm;  Vec              coordinates;  const PetscReal *quadPoints;  PetscScalar     *vertices = NULL;  PetscInt         numQuadPoints, csize, dim, d, q;  PetscErrorCode   ierr;  PetscFunctionBegin;  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);  ierr = DMDAGetInfo(dm, &dim, 0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);  ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);  ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);  ierr = DMDAVecGetClosure(cdm, NULL, coordinates, cell, &csize, &vertices);CHKERRQ(ierr);  for (d = 0; d < dim; ++d) v0[d] = PetscRealPart(vertices[d]);  switch (dim) {  case 2:    ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, NULL);CHKERRQ(ierr);    for (q = 0; q < numQuadPoints; ++q) {      ierr = DMDAComputeCellGeometry_2D(dm, vertices, &quadPoints[q*dim], J, invJ, detJ);CHKERRQ(ierr);    }    break;  default:    SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Dimension %d not supported", dim);  }  ierr = DMDAVecRestoreClosure(cdm, NULL, coordinates, cell, &csize, &vertices);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode SetCoordinates3d(DM da){  PetscErrorCode ierr;  PetscInt       i,j,mstart,m,nstart,n,pstart,p,k;  Vec            gc,global;  DMDACoor3d     ***coors;  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 = DMGetCoordinatesLocal(da,&gc);CHKERRQ(ierr);  ierr = DMDAVecGetArray(cda,gc,&coors);CHKERRQ(ierr);  ierr = DMDAGetCorners(cda,&mstart,&nstart,&pstart,&m,&n,&p);CHKERRQ(ierr);  for (i=mstart; i<mstart+m; i++) {    for (j=nstart; j<nstart+n; j++) {      for (k=pstart; k<pstart+p; k++) {        if (i % 2) {          coors[k][j][i].x = coors[k][j][i-1].x + .1*(coors[k][j][i+1].x - coors[k][j][i-1].x);        }        if (j % 2) {          coors[k][j][i].y = coors[k][j-1][i].y + .3*(coors[k][j+1][i].y - coors[k][j-1][i].y);        }        if (k % 2) {          coors[k][j][i].z = coors[k-1][j][i].z + .4*(coors[k+1][j][i].z - coors[k-1][j][i].z);        }      }    }  }  ierr = DMDAVecRestoreArray(cda,gc,&coors);CHKERRQ(ierr);  ierr = DMGetCoordinates(da,&global);CHKERRQ(ierr);  ierr = DMLocalToGlobalBegin(cda,gc,INSERT_VALUES,global);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(cda,gc,INSERT_VALUES,global);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode SetCoordinates1d(DM da){  PetscErrorCode ierr;  PetscInt       i,start,m;  Vec            gc,global;  PetscScalar    *coors;  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 = DMGetCoordinatesLocal(da,&gc);CHKERRQ(ierr);  ierr = DMDAVecGetArray(cda,gc,&coors);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] = coors[i-1] + .1*(coors[i+1] - coors[i-1]);    }  }  ierr = DMDAVecRestoreArray(cda,gc,&coors);CHKERRQ(ierr);  ierr = DMGetCoordinates(da,&global);CHKERRQ(ierr);  ierr = DMLocalToGlobalBegin(cda,gc,INSERT_VALUES,global);CHKERRQ(ierr);  ierr = DMLocalToGlobalEnd(cda,gc,INSERT_VALUES,global);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 ini_bou(Vec X,AppCtx* user){  PetscErrorCode ierr;  DM             cda;  DMDACoor2d     **coors;  PetscScalar    **p;  Vec            gc;  PetscInt       M,N,I,J;  PetscMPIInt    rank;  PetscFunctionBeginUser;  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);  ierr = DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);  user->dx = (user->xmax - user->xmin)/(M-1); user->dy = (user->ymax - user->ymin)/(N-1);  ierr = DMGetCoordinateDM(user->da,&cda);CHKERRQ(ierr);  ierr = DMGetCoordinates(user->da,&gc);CHKERRQ(ierr);  ierr = DMDAVecGetArrayRead(cda,gc,&coors);CHKERRQ(ierr);  ierr = DMDAVecGetArray(user->da,X,&p);CHKERRQ(ierr);  /* Point mass at (mux,muy) */  ierr = PetscPrintf(PETSC_COMM_WORLD,"Original user->mux = %f, user->muy = %f/n",user->mux,user->muy);CHKERRQ(ierr);  ierr = DMDAGetLogicalCoordinate(user->da,user->mux,user->muy,0.0,&I,&J,NULL,&user->mux,&user->muy,NULL);CHKERRQ(ierr);  user->PM_min = user->Pmax*PetscSinScalar(user->mux);  ierr = PetscPrintf(PETSC_COMM_WORLD,"Corrected user->mux = %f, user->muy = %f user->PM_min = %f,user->dx = %f/n",user->mux,user->muy,user->PM_min,user->dx);CHKERRQ(ierr);  if (I > -1 && J > -1) {    p[J][I] = 1.0;  }  ierr = DMDAVecRestoreArrayRead(cda,gc,&coors);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(user->da,X,&p);CHKERRQ(ierr);  PetscFunctionReturn(0);}

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

/*  FormExactSolution2 - Forms initial approximation.  Input Parameters:  da - The DM  user - user-defined application context  Output Parameter:  X - vector */PetscErrorCode FormExactSolution2(DM da, AppCtx *user, Vec U){  DM             coordDA;  Vec            coordinates;  DMDACoor2d   **coords;  PetscScalar  **u;  PetscReal      x, y;  PetscInt       xs, ys, xm, ym, i, j;  PetscErrorCode ierr;  PetscFunctionBeginUser;  ierr = DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);CHKERRQ(ierr);  ierr = DMGetCoordinateDM(da, &coordDA);CHKERRQ(ierr);  ierr = DMGetCoordinates(da, &coordinates);CHKERRQ(ierr);  ierr = DMDAVecGetArray(coordDA, coordinates, &coords);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da, U, &u);CHKERRQ(ierr);  for (j = ys; j < ys+ym; ++j) {    for (i = xs; i < xs+xm; ++i) {      x = PetscRealPart(coords[j][i].x);      y = PetscRealPart(coords[j][i].y);      u[j][i] = PetscSinReal(PETSC_PI*x)*PetscSinReal(PETSC_PI*y);    }  }  ierr = DMDAVecRestoreArray(da, U, &u);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(coordDA, coordinates, &coords);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode FormPsiAndExactSoln(DM da) {  ObsCtx         *user;  PetscErrorCode ierr;  DMDALocalInfo  info;  PetscInt       i,j;  DM             coordDA;  Vec            coordinates;  DMDACoor2d     **coords;  PetscReal      **psi, **uexact, r;  const PetscReal afree = 0.69797, A = 0.68026, B = 0.47152;  PetscFunctionBeginUser;  ierr = DMGetApplicationContext(da,&user);CHKERRQ(ierr);  ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);  ierr = DMGetCoordinateDM(da, &coordDA);CHKERRQ(ierr);  ierr = DMGetCoordinates(da, &coordinates);CHKERRQ(ierr);  ierr = DMDAVecGetArray(coordDA, coordinates, &coords);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da, user->psi, &psi);CHKERRQ(ierr);  ierr = DMDAVecGetArray(da, user->uexact, &uexact);CHKERRQ(ierr);  for (j=info.ys; j<info.ys+info.ym; j++) {    for (i=info.xs; i<info.xs+info.xm; i++) {      r = PetscSqrtReal(pow(coords[j][i].x,2) + pow(coords[j][i].y,2));      if (r <= 1.0) psi[j][i] = PetscSqrtReal(1.0 - r * r);      else psi[j][i] = -1.0;      if (r <= afree) uexact[j][i] = psi[j][i];  /* on the obstacle */      else uexact[j][i] = - A * PetscLogReal(r) + B;   /* solves the laplace eqn */    }  }  ierr = DMDAVecRestoreArray(da, user->psi, &psi);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(da, user->uexact, &uexact);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(coordDA, coordinates, &coords);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);}

PetscErrorCode ini_bou(Vec X,AppCtx* user){  PetscErrorCode ierr;  DM             cda;  DMDACoor2d     **coors;  PetscScalar    **p;  Vec            gc;  PetscInt       i,j;  PetscInt       xs,ys,xm,ym,M,N;  PetscScalar    xi,yi;  PetscScalar    sigmax=user->sigmax,sigmay=user->sigmay;  PetscScalar    rho   =user->rho;  PetscScalar    mux   =user->mux,muy=user->muy;  PetscMPIInt    rank;  PetscScalar    sum;  PetscFunctionBeginUser;  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);  ierr = DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);  user->dx = (user->xmax - user->xmin)/(M-1); user->dy = (user->ymax - user->ymin)/(N-1);  ierr = DMGetCoordinateDM(user->da,&cda);CHKERRQ(ierr);  ierr = DMGetCoordinates(user->da,&gc);CHKERRQ(ierr);  ierr = DMDAVecGetArray(cda,gc,&coors);CHKERRQ(ierr);  ierr = DMDAVecGetArray(user->da,X,&p);CHKERRQ(ierr);  ierr = DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);  /* mux and muy need to be grid points in the x and y-direction otherwise the solution goes unstable     muy is set by choosing the y domain, no. of grid points along y-direction so that muy is a grid point     in the y-direction. We only modify mux here  */  mux = user->mux = coors[0][M/2+10].x; /* For -pi < x < pi, this should be some angle between 0 and pi/2 */  if (user->nonoiseinitial) {    for (i=xs; i < xs+xm; i++) {      for (j=ys; j < ys+ym; j++) {        xi = coors[j][i].x; yi = coors[j][i].y;        if ((xi == mux) && (yi == muy)) {          p[j][i] = 1.0;        }      }    }  } else {    /* Change PM_min accordingly */    user->PM_min = user->Pmax*sin(mux);    for (i=xs; i < xs+xm; i++) {      for (j=ys; j < ys+ym; j++) {        xi = coors[j][i].x; yi = coors[j][i].y;        p[j][i] = (0.5/(PETSC_PI*sigmax*sigmay*PetscSqrtScalar(1.0-rho*rho)))*PetscExpScalar(-0.5/(1-rho*rho)*(PetscPowScalar((xi-mux)/sigmax,2) + PetscPowScalar((yi-muy)/sigmay,2) - 2*rho*(xi-mux)*(yi-muy)/(sigmax*sigmay)));      }    }  }  ierr = DMDAVecRestoreArray(cda,gc,&coors);CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(user->da,X,&p);CHKERRQ(ierr);  ierr = VecSum(X,&sum);CHKERRQ(ierr);  ierr = VecScale(X,1.0/sum);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode DAApplyTrilinearMapping(DM da){  PetscErrorCode ierr;  PetscInt       i,j,k;  PetscInt       sx,nx,sy,ny,sz,nz;  Vec            Gcoords;  DMDACoor3d     ***XX;  PetscScalar    xx,yy,zz;  DM             cda;  PetscFunctionBeginUser;  ierr = DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0 );CHKERRQ(ierr);  ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr);  ierr = DMGetCoordinates(da,&Gcoords);CHKERRQ(ierr);  ierr = DMDAVecGetArray(cda,Gcoords,&XX);CHKERRQ(ierr);  ierr = DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz);CHKERRQ(ierr);  for (i=sx; i<sx+nx; i++) {    for (j=sy; j<sy+ny; j++ ) {      for (k=sz; k<sz+nz; k++ ) {        PetscScalar Ni[8];        PetscScalar xi   = XX[k][j][i].x;        PetscScalar eta  = XX[k][j][i].y;        PetscScalar zeta = XX[k][j][i].z;        PetscScalar xn[] = {0.0,2.0,0.2,3.5,   0.0,2.1,0.23,3.125  };        PetscScalar yn[] = {-1.3,0.0,2.0,4.0,  -1.45,-0.1,2.24,3.79  };        PetscScalar zn[] = {0.0,0.3,-0.1,0.123,  0.956,1.32,1.12,0.798  };        PetscInt    p;        Ni[0] = 0.125*(1.0-xi)*(1.0-eta)*(1.0-zeta);        Ni[1] = 0.125*(1.0+xi)*(1.0-eta)*(1.0-zeta);        Ni[2] = 0.125*(1.0-xi)*(1.0+eta)*(1.0-zeta);        Ni[3] = 0.125*(1.0+xi)*(1.0+eta)*(1.0-zeta);        Ni[4] = 0.125*(1.0-xi)*(1.0-eta)*(1.0+zeta);        Ni[5] = 0.125*(1.0+xi)*(1.0-eta)*(1.0+zeta);        Ni[6] = 0.125*(1.0-xi)*(1.0+eta)*(1.0+zeta);        Ni[7] = 0.125*(1.0+xi)*(1.0+eta)*(1.0+zeta);        xx = yy = zz = 0.0;        for (p=0; p<8; p++ ) {          xx += Ni[p]*xn[p];          yy += Ni[p]*yn[p];          zz += Ni[p]*zn[p];        }        XX[k][j][i].x = xx;        XX[k][j][i].y = yy;        XX[k][j][i].z = zz;      }    }  }  ierr = DMDAVecRestoreArray(cda,Gcoords,&XX);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);}

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

PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat *J,Mat *Jpre,MatStructure *flg,void *ctx){  PetscErrorCode ierr;  AppCtx         *user=(AppCtx*)ctx;  DM             cda;  DMDACoor2d     **coors;  PetscInt       i,j;  PetscInt       xs,ys,xm,ym,M,N;  Vec            gc;  PetscScalar    val[5],xi,yi;  MatStencil     row,col[5];  PetscScalar    c1,c3,c5,c1pos,c1neg,c3pos,c3neg;  PetscFunctionBeginUser;  *flg = SAME_NONZERO_PATTERN;  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 = DMGetCoordinatesLocal(user->da,&gc);CHKERRQ(ierr);  ierr = DMDAVecGetArray(cda,gc,&coors);CHKERRQ(ierr);  for (i=xs; i < xs+xm; i++) {    for (j=ys; j < ys+ym; j++) {      PetscInt nc = 0;      xi = coors[j][i].x; yi = coors[j][i].y;      row.i = i; row.j = j;      c1        = (yi-user->ws)/user->dx;      c1pos    = PetscMax(c1,0);      c1neg    = PetscMin(c1,0);      c3        = (user->ws/(2.0*user->H))*(user->PM_min - user->Pmax*sin(xi))/user->dy;      c3pos    = PetscMax(c3,0);      c3neg    = PetscMin(c3,0);      c5        = (PetscPowScalar((user->lambda*user->ws)/(2*user->H),2)*user->q*(1.0-PetscExpScalar(-t/user->lambda)))/(user->dy*user->dy);      col[nc].i = i-1; col[nc].j = j;   val[nc++] = c1pos;      col[nc].i = i+1; col[nc].j = j;   val[nc++] = -c1neg;      col[nc].i = i;   col[nc].j = j-1; val[nc++] = c3pos + c5;      col[nc].i = i;   col[nc].j = j+1; val[nc++] = -c3neg + c5;      col[nc].i = i;   col[nc].j = j;   val[nc++] = -c1pos + c1neg -c3pos + c3neg -2*c5 -a;      ierr = MatSetValuesStencil(*Jpre,1,&row,nc,col,val,INSERT_VALUES);CHKERRQ(ierr);    }  }  ierr = DMDAVecRestoreArray(cda,gc,&coors);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);  }  PetscFunctionReturn(0);}

/*@C   DMDAGetLogicalCoordinate - Returns a the i,j,k logical coordinate for the closest mesh point to a x,y,z point in the coordinates of the DMDA   Collective on DMDA   Input Parameters:+  da - the distributed array-  x,y,z - the physical coordinates   Output Parameters:+   II, JJ, KK - the logical coordinate (-1 on processes that do not contain that point)-   X, Y, Z, - (optional) the coordinates of the located grid point   Level: advanced   Notes:   All processors that share the DMDA must call this with the same coordinate value.keywords: distributed array, get, processor subset@*/PetscErrorCode  DMDAGetLogicalCoordinate(DM da,PetscScalar x,PetscScalar y,PetscScalar z,PetscInt *II,PetscInt *JJ,PetscInt *KK,PetscScalar *X,PetscScalar *Y,PetscScalar *Z){  DM_DA          *dd = (DM_DA*)da->data;  PetscErrorCode ierr;  Vec            coors;  DM             dacoors;  DMDACoor2d     **c;  PetscInt       i,j,xs,xm,ys,ym;  PetscReal      d,D = PETSC_MAX_REAL,Dv;  PetscMPIInt    rank,root;  PetscFunctionBegin;  if (dd->dim == 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot get point from 1d DMDA");  if (dd->dim == 3) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot get point from 3d DMDA");  *II = -1;  *JJ = -1;  ierr = DMGetCoordinateDM(da,&dacoors);CHKERRQ(ierr);  ierr = DMDAGetCorners(dacoors,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);  ierr = DMGetCoordinates(da,&coors);CHKERRQ(ierr);  ierr = DMDAVecGetArray(dacoors,coors,&c);CHKERRQ(ierr);  for (j=ys; j<ys+ym; j++) {    for (i=xs; i<xs+xm; i++) {      d = PetscSqrtReal(PetscRealPart( (c[j][i].x - x)*(c[j][i].x - x) + (c[j][i].y - y)*(c[j][i].y - y) ));      if (d < D) {        D   = d;        *II = i;        *JJ = j;      }    }  }  ierr = MPI_Allreduce(&D,&Dv,1,MPIU_REAL,MPI_MIN,PetscObjectComm((PetscObject)da));CHKERRQ(ierr);  if (D != Dv) {    *II  = -1;    *JJ  = -1;    rank = 0;  } else {    *X = c[*JJ][*II].x;    *Y = c[*JJ][*II].y;    ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)da),&rank);CHKERRQ(ierr);    rank++;  }  ierr = MPI_Allreduce(&rank,&root,1,MPI_INT,MPI_SUM,PetscObjectComm((PetscObject)da));CHKERRQ(ierr);  root--;  ierr = MPI_Bcast(X,1,MPIU_SCALAR,root,PetscObjectComm((PetscObject)da));CHKERRQ(ierr);  ierr = MPI_Bcast(Y,1,MPIU_SCALAR,root,PetscObjectComm((PetscObject)da));CHKERRQ(ierr);  ierr = DMDAVecRestoreArray(dacoors,coors,&c);CHKERRQ(ierr);  PetscFunctionReturn(0);}

/* FormFunctionLocalMMS2 - Evaluates nonlinear function, F(x) on local process patch */PetscErrorCode FormFunctionLocalMMS2(DMDALocalInfo *info,PetscScalar **vx,PetscScalar **f,AppCtx *user){  PetscErrorCode ierr;  PetscInt       i,j;  PetscReal      lambda,hx,hy,hxdhy,hydhx;  PetscScalar    u,ue,uw,un,us,uxx,uyy;  PetscReal      x,y;  DM             coordDA;  Vec            coordinates;  DMDACoor2d   **coords;  PetscFunctionBeginUser;  lambda = user->param;  hx     = 1.0/(PetscReal)(info->mx-1);  hy     = 1.0/(PetscReal)(info->my-1);  hxdhy  = hx/hy;  hydhx  = hy/hx;  /* Extract coordinates */  ierr = DMGetCoordinateDM(info->da, &coordDA);CHKERRQ(ierr);  ierr = DMGetCoordinates(info->da, &coordinates);CHKERRQ(ierr);  ierr = DMDAVecGetArray(coordDA, coordinates, &coords);CHKERRQ(ierr);  /* Compute function over the locally owned part of the grid */  for (j=info->ys; j<info->ys+info->ym; j++) {    for (i=info->xs; i<info->xs+info->xm; i++) {      if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {        f[j][i] = 2.0*(hydhx+hxdhy)*vx[j][i];      } else {        x  = PetscRealPart(coords[j][i].x);        y  = PetscRealPart(coords[j][i].y);        u  = vx[j][i];        uw = vx[j][i-1];        ue = vx[j][i+1];        un = vx[j-1][i];        us = vx[j+1][i];        if (i-1 == 0) uw = 0.;        if (i+1 == info->mx-1) ue = 0.;        if (j-1 == 0) un = 0.;        if (j+1 == info->my-1) us = 0.;        uxx     = (2.0*u - uw - ue)*hydhx;        uyy     = (2.0*u - un - us)*hxdhy;        f[j][i] = uxx + uyy - hx*hy*(lambda*PetscExpScalar(u) + 2*PetscSqr(PETSC_PI)*PetscSinReal(PETSC_PI*x)*PetscSinReal(PETSC_PI*y) - lambda*exp(PetscSinReal(PETSC_PI*x)*PetscSinReal(PETSC_PI*y)));      }    }  }  ierr = DMDAVecRestoreArray(coordDA, coordinates, &coords);CHKERRQ(ierr);  ierr = PetscLogFlops(11.0*info->ym*info->xm);CHKERRQ(ierr);  PetscFunctionReturn(0);}

PetscErrorCode ini_bou(Vec X,AppCtx* user){    PetscErrorCode ierr;    DM             cda;    DMDACoor2d     **coors;    PetscScalar    **p;    Vec            gc;    PetscInt       i,j;    PetscInt       xs,ys,xm,ym,M,N;    PetscScalar    xi,yi;    PetscScalar    sigmax=user->sigmax,sigmay=user->sigmay;    PetscScalar    rho   =user->rho;    PetscScalar    mux   =user->mux,muy=user->muy;    PetscMPIInt    rank;    PetscFunctionBeginUser;    ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);    CHKERRQ(ierr);    ierr = DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);    user->dx = (user->xmax - user->xmin)/(M-1);    user->dy = (user->ymax - user->ymin)/(N-1);    ierr = DMGetCoordinateDM(user->da,&cda);    CHKERRQ(ierr);    ierr = DMGetCoordinates(user->da,&gc);    CHKERRQ(ierr);    ierr = DMDAVecGetArray(cda,gc,&coors);    CHKERRQ(ierr);    ierr = DMDAVecGetArray(user->da,X,&p);    CHKERRQ(ierr);    ierr = DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);    CHKERRQ(ierr);    for (i=xs; i < xs+xm; i++) {        for (j=ys; j < ys+ym; j++) {            xi = coors[j][i].x;            yi = coors[j][i].y;            if (i == 0 || j == 0 || i == M-1 || j == N-1) p[j][i] = 0.0;            else p[j][i] = (0.5/(PETSC_PI*sigmax*sigmay*PetscSqrtScalar(1.0-rho*rho)))*PetscExpScalar(-0.5/(1-rho*rho)*(PetscPowScalar((xi-mux)/sigmax,2) + PetscPowScalar((yi-muy)/sigmay,2) - 2*rho*(xi-mux)*(yi-muy)/(sigmax*sigmay)));        }    }    /*  p[N/2+N%2][M/2+M%2] = 1/(user->dx*user->dy); */    ierr = DMDAVecRestoreArray(cda,gc,&coors);    CHKERRQ(ierr);    ierr = DMDAVecRestoreArray(user->da,X,&p);    CHKERRQ(ierr);    PetscFunctionReturn(0);}