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: }