Actual source code: ex29.c

  2: static char help[] = "Tests DA wirebasket interpolation.\n\n";

 4:  #include petscda.h


 11: int main(int argc,char **argv)
 12: {
 14:   DA             da;
 15:   Mat            Aglobal,P;

 17:   PetscInitialize(&argc,&argv,(char*)0,help);

 19:   DACreate3d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX,3,5,5,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,PETSC_NULL,&da);
 20:   DAGetMatrix(da,MATAIJ,&Aglobal);
 21:   ComputeMatrix(da,Aglobal);
 22:   DAGetWireBasketInterpolation(da,Aglobal,MAT_INITIAL_MATRIX,&P);
 23:   MatView(P,PETSC_VIEWER_STDOUT_WORLD);

 25:   MatDestroy(P);
 26:   MatDestroy(Aglobal);
 27:   DADestroy(da);
 28:   PetscFinalize();
 29:   return 0;
 30: }
 31: 
 34: PetscErrorCode ComputeMatrix(DA da,Mat B)
 35: {
 37:   PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,cnt;
 38:   PetscScalar    v[7],Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy,val;
 39:   MatStencil     row,col[7];

 41:   DAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0);
 42:   Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
 43:   HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
 44:   DAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
 45: 
 46: #define Endpoint(a,b) (a == 0 || a == (b-1))

 48:   for (k=zs; k<zs+zm; k++){
 49:     for (j=ys; j<ys+ym; j++){
 50:       for(i=xs; i<xs+xm; i++){
 51:         row.i = i; row.j = j; row.k = k;
 52:         cnt   = 0;
 53:         val   = 0;
 54:         if (k != 0) { val += v[cnt] = -HxHydHz;col[cnt].i = i; col[cnt].j = j; col[cnt].k = k-1;cnt++;}
 55:         if (j != 0) { val += v[cnt] = -HxHzdHy;col[cnt].i = i; col[cnt].j = j-1; col[cnt].k = k;cnt++;}
 56:         if (i != 0) { val += v[cnt] = -HyHzdHx;col[cnt].i = i-1; col[cnt].j = j; col[cnt].k = k;cnt++;}
 57:         if (i != mx-1) {val += v[cnt] = -HyHzdHx;col[cnt].i = i+1; col[cnt].j = j; col[cnt].k = k;cnt++;}
 58:         if (j != my-1) {val += v[cnt] = -HxHzdHy;col[cnt].i = i; col[cnt].j = j+1; col[cnt].k = k;cnt++;}
 59:         if (k != mz-1) {val += v[cnt] = -HxHydHz;col[cnt].i = i; col[cnt].j = j; col[cnt].k = k+1;cnt++;}
 60:         v[cnt] = -val;col[cnt].i = row.i; col[cnt].j = row.j; col[cnt].k = row.k;cnt++;
 61:         MatSetValuesStencil(B,1,&row,cnt,col,v,INSERT_VALUES);
 62:       }
 63:     }
 64:   }
 65:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 66:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
 67:   return 0;
 68: }