Actual source code: ex29.c
2: static char help[] = "Tests DA wirebasket interpolation.\n\n";
4: #include petscda.h
5: #include petscsys.h
12: int main(int argc,char **argv)
13: {
15: DA da;
16: Mat Aglobal,P;
18: PetscInitialize(&argc,&argv,(char*)0,help);
20: 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);
21: DAGetMatrix(da,MATAIJ,&Aglobal);
22: ComputeMatrix(da,Aglobal);
23: DAGetWireBasketInterpolation(da,Aglobal,MAT_INITIAL_MATRIX,&P);
24: MatView(P,PETSC_VIEWER_STDOUT_WORLD);
26: MatDestroy(P);
27: MatDestroy(Aglobal);
28: DADestroy(da);
29: PetscFinalize();
30: return 0;
31: }
32:
35: PetscErrorCode ComputeMatrix(DA da,Mat B)
36: {
38: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,cnt;
39: PetscScalar v[7],Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy,val;
40: MatStencil row,col[7];
42: DAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0);
43: Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
44: HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
45: DAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
46:
47: #define Endpoint(a,b) (a == 0 || a == (b-1))
49: for (k=zs; k<zs+zm; k++){
50: for (j=ys; j<ys+ym; j++){
51: for(i=xs; i<xs+xm; i++){
52: row.i = i; row.j = j; row.k = k;
53: cnt = 0;
54: val = 0;
55: if (k != 0) { val += v[cnt] = -HxHydHz;col[cnt].i = i; col[cnt].j = j; col[cnt].k = k-1;cnt++;}
56: if (j != 0) { val += v[cnt] = -HxHzdHy;col[cnt].i = i; col[cnt].j = j-1; col[cnt].k = k;cnt++;}
57: if (i != 0) { val += v[cnt] = -HyHzdHx;col[cnt].i = i-1; col[cnt].j = j; col[cnt].k = k;cnt++;}
58: if (i != mx-1) {val += v[cnt] = -HyHzdHx;col[cnt].i = i+1; col[cnt].j = j; col[cnt].k = k;cnt++;}
59: if (j != my-1) {val += v[cnt] = -HxHzdHy;col[cnt].i = i; col[cnt].j = j+1; col[cnt].k = k;cnt++;}
60: if (k != mz-1) {val += v[cnt] = -HxHydHz;col[cnt].i = i; col[cnt].j = j; col[cnt].k = k+1;cnt++;}
61: v[cnt] = -val;col[cnt].i = row.i; col[cnt].j = row.j; col[cnt].k = row.k;cnt++;
62: MatSetValuesStencil(B,1,&row,cnt,col,v,INSERT_VALUES);
63: }
64: }
65: }
66: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
67: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
68: return 0;
69: }