[petsc-users] changing DAs matrix and coefficients ordering

Jed Brown jed at 59A2.org
Tue Feb 2 11:10:50 CST 2010


On Tue, 2 Feb 2010 17:36:34 +0100 (CET), "ilyas yilmaz" <ilyasy at chalmers.se> wrote:
> I have actually tried it before, but I failed. I guess I am missing
> something. Could you please explain me how to do this on a simple DAs
> example?

Suppose you have an M by N region that you want to index as x[i][j] with
i in [0..M) and j in [0..N).

  DACreate2d(comm,DA_NOPERIODIC,DA_STENCIL_BOX,N,M,n,m,dof,s,ly,lx,&da);

then something like

  PetscErrorCode FormFunction(DA da,Vec X,Vec F,void *ctx) {
    Node **x,**f;
    PetscInt i,j,xs,ys,xm,ym,mx,my;
    ...
    DAGetInfo(da,0, &my,&mx,0, 0,0,0, 0,0,0,0);
    hx = Lx / mx;
    hy = Ly / my;
    DAGetCorners(da,&ys,&xs,0,&ym,&xm,0);
    DAVecGetArray(da,X,&x);
    DAVecGetArray(da,F,&f);
    for (i=xs; i<xs+xm; i++) {
      for (j=ys; j<ys+ym; j++) {
        f[i][j].u = -x[i-1][j].u -x[i][j-1].u + 4*x[i][j].u - x[i][j+1].u - x[i+1][j].u;
        ...


If you're still confused and want a working example, PETSc-dev has

  http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-dev/src/snes/examples/tutorials/ex48.c.html

but it is overly complex to demonstrate this point since it manages both
a 2D and 3D DA with compatible layouts.

Jed


More information about the petsc-users mailing list