petsc 结构化数据使用
前文讲到petsc两类基本数据对象（mat, vec)。而petsc的这两类数据对象默认都是按一维存储和创建的。实际pde问题，往往是2维或者3维的。故petsc内建一些映射函数，将默认的一维数据映射成逻辑上的多维数组。本文主要介绍向规则的2D数组映射的函数。
petsc中逻辑上的多维数组称为 "DA" ( distributed array), 创建一个2D DA：
DACreate2d(MPI_Comm comm, DAPeriodicType wrap, DAStencilType stencil_type, PetscInt M, PetscInt N, PetscINt m, PetscInt n, PetscInt dof, PetscInt s, PetscInt *lx, PetscInt *ly, DA *da)
/* comm MPI通信容器
wrap 周期边界
stencil_type 差分模式（box, star)
M ， N x, y方向整体空间节点数
m, n ， x, y方向cpu节点数
dof 空间节点自由度数
s 差分模式的厚度(即，how many 层数 空间节点）
*lx, *ly cpu节点所管辖的部分节点数组
*da, 所创建DA对象 * 获取该DA对象的全局及局部信息（即分配到cpu节点上的单元信息）采用一 DAGetLocalInfo (DA da, DALocalInfo *info) ，其中
typedef struct {
PetscInt dim, dof, sw ;
PetscInt mx, my, mz; /* global number of grid points in x, y, z */
PetscInt xs, ys, zs; /* indices to the start of grid points in current cpu, excluiding ghost points(gp) */
PetscInt xm, ym, zm;/* local number of grid points in current cpu, no GP */
PetscInt gxs, gys, gzs; /* similar to xs, but including GP */
PetscInt gxm, gym, gzm; /* similar to xm, but including GP */
DAPeriodicType pt;
DAStencilType st;
DA da;
} DALocalInfo ;DA对象相当于一个逻辑层，用于管理其所有空间节点以及该通信器里的所有cpu节点。与之对应的物理层，就是DAVec 或者 DAMat。 将数据通过逻辑层管理的好处是直观。 比如，将一个计算域离散成n * n 单元 ，采用DA表述就是我们头脑中的 n * n 的样子，而在物理层仍然是一维线，不直观。
但是，数据实际上仍然是存在物理层的。因此对数据操作还是要再回到VEC， MAT。 （抽象出来的DA一方面是更直观了，但是又不能直接使用 ）
正确地讲，DA就是规则网格信息。DAVectGetArray函数用于建立列向量，且该列向量元素个数等于DA网格节点数×自由度数。DAGetMatrix用于建立矩阵A，且矩阵宽度等于DA两个维度之积。
DACreateGlobalVector(DA da, Vec *global_vec)
DACreateLocalVector(DA da, Vec *local_vec)
/* 分别声明了一个全局向量和一个本地向量(包含GP），均基于DA对象。*/
DAVecGetArray (DA da, Vec vec, void* array) ; /* 赋值需要对物理层的VEC进行 */
DAVecRestoreArray(DA da, Vec vec, void* array) ; /* 赋值完成，再回到逻辑层DA ×/
举例：全局（整体）解向量 u(x,y ) =0 , v(x, y) = 1。首先创建DA，（造了空房子）；然后通过DA创建GlobalVector，（将空房子与一个GLOBALVECT关联起来 ）； 然后获取一个指向该DA（已包含将要赋值的GLOBALVECT）的指针，再通过该指针完成对元素的访问/修改。
DA da;
DALocalInfo info;
PetscInt i,j;
Vec solv;
typedef struct {
PEtscScalar u, v;
} Field ;
Field **g;
DACreate2d(...);
DACreateGlobalVector(da, &solv);
DAGetLocalInfo(da, &info);
DAVecGetArray(da, solv, &g);
for(i=info.xs; i < info.xs + info.xm; i++ ){
for(j=info.ys; j<info.ys + info.ym; j++){
g[j][i].u = 0;
g[j][i].v = 1;
}
}
DAVecRestoreArray(da, solv, &g);另外，更新GP采用：
DAGlobalToLocalBegin(DA da, Vec global, InsertMode mode, Vec local);
DAGlobalToLocalend(DA da, Vec global, InsertMode mode, Vec local); 类似的，DA 到MAT 的映射：
DAGetMatrix(DA da, MatType mtype, Mat *mat) /* create Matrix based on DA array */
MatSetValuesStencil(Mat mat, PetscInt m, const MatStencil row[], PetscINt n, const MatStencil col[], const PetscScalar val[], InsertMode mode)
/*
m, n :: number of rows and colums;
row[], col[] :: grid indices(associated with stencil structure(pattern)); and components indices;
val :: values to insert
petsc默认按行划分，col[]实际上就存当前行中不为0的元素。
*/
/* 对MAT赋值完之后，类似VEC，需要再返映射到“逻辑层”*/
MatAssemblyBegin(Mat mat, MatAssemblyType type);
MatAssemblyEnd(Mat mat, MatAssemblyType type);
举例：一维laplace差分离散
DA da; DALocalInfo info; Mat M; PetscInt i,N=25, num_cols; MatStructure row, col[3]; PetscScalar val[3], h=1.00/(N-1); DACreate1d(...); DAGetLocalInfo(da, &info); DAGetMatrix(da, MATMPIAIJ, &M); for(i=info.xs; i < info.xs + info.xm; i++ ){ row.i = i ; /* grid space indices */ row.c = 0 ; /* first component */ if(i ==0 || i = N -1 ) /* B.C. */ { col[0].i = i ; col[0].c = 0 ; /* only one component */ val[0] = 1.00 ; num_cols = 1; } else { col[0].i = i -1;/* left point */ col[0].c = 0; val[0] = 1/h/h; col[1].i = i; col[1].c = 0; val[1] = -2/h/h; col[2].i = i+1; col[2].c = 0; val[2] = 1/h/h; num_cols = 3; } MatSetValuesStencil(M,1,&row, num_cols, col, val, INSERT_VALUES); }MatAssemblyBegin(Mat mat, MatAssemblyType type); MatAssemblyEnd(Mat mat, MatAssemblyType type);
总结： 规则计算域的数据采用DA管理较直观 （映射成DA）；但DA对象不能直接访问/赋值VEC/MAT数据。故需将DA对象返映射回VECT/MAT。然后对于VECT，使用函数DAVecGetArray返回指向该VECT的指针，并由该指针访问/赋值VECT的元素；对于MAT，使用函数MatSetValuesStencil实现访问。
DA管理规则计算域，关联此计算域的方程组A * x = b 的矩阵A，列向量b 都可由DA获取。
实际上，DA对象很绕圈子很麻烦。但这正是petsc的设计原则之一：petsc函数对用户数据不可见，所以只好通过指针访问。