xref: /petsc/src/dm/impls/da/dagetarray.c (revision e611a964e9853b74d61a56642fe9d06a6e51780f) 
1 
2 #include <petsc/private/dmdaimpl.h>    /*I   "petscdmda.h"   I*/
3 
4 #undef __FUNCT__
5 #define __FUNCT__ "DMDAVecGetArray"
6 /*@C
7    DMDAVecGetArray - Returns a multiple dimension array that shares data with
8       the underlying vector and is indexed using the global dimensions.
9 
10    Logically collective on Vec
11 
12    Input Parameter:
13 +  da - the distributed array
14 -  vec - the vector, either a vector the same size as one obtained with DMCreateGlobalVector() or DMCreateLocalVector()
15 
16    Output Parameter:
17 .  array - the array
18 
19    Notes:
20     Call DMDAVecRestoreArray() once you have finished accessing the vector entries.
21 
22     In C, the indexing is "backwards" from what expects: array[k][j][i] NOT array[i][j][k]!
23 
24     If vec is a local vector (obtained with DMCreateLocalVector() etc) then the ghost point locations are accessible. If it is
25     a global vector then the ghost points are not accessible. Of course with the local vector you will have had to do the
26 
27     appropriate DMGlobalToLocalBegin() and DMGlobalToLocalEnd() to have correct values in the ghost locations.
28 
29   Fortran Notes: From Fortran use DMDAVecGetArrayF90() and pass for the array type PetscScalar,pointer :: array(:,...,:) of the appropriate
30        dimension. For a DMDA created with a dof of 1 use the dimension of the DMDA, for a DMDA created with a dof greater than 1 use one more than the
31        dimension of the DMDA. The order of the indices is array(xs:xs+xm-1,ys:ys+ym-1,zs:zs+zm-1) (when dof is 1) otherwise
32        array(0:dof-1,xs:xs+xm-1,ys:ys+ym-1,zs:zs+zm-1) where the values are obtained from
33        DMDAGetCorners() for a global array or DMDAGetGhostCorners() for a local array. Include petsc/finclude/petscdmda.h90 to access this routine.
34 
35   Due to bugs in the compiler DMDAVecGetArrayF90() does not work with gfortran versions before 4.5
36 
37   Level: intermediate
38 
39 .keywords: distributed array, get, corners, nodes, local indices, coordinates
40 
41 .seealso: DMDAGetGhostCorners(), DMDAGetCorners(), VecGetArray(), VecRestoreArray(), DMDAVecRestoreArray(), DMDAVecRestoreArrayDOF()
42           DMDAVecGetArrayDOF()
43 @*/
44 PetscErrorCode  DMDAVecGetArray(DM da,Vec vec,void *array)
45 {
46   PetscErrorCode ierr;
47   PetscInt       xs,ys,zs,xm,ym,zm,gxs,gys,gzs,gxm,gym,gzm,N,dim,dof;
48 
49   PetscFunctionBegin;
50   PetscValidHeaderSpecific(da, DM_CLASSID, 1);
51   PetscValidHeaderSpecific(vec, VEC_CLASSID, 2);
52   PetscValidPointer(array, 3);
53   if (da->defaultSection) {
54     ierr = VecGetArray(vec,(PetscScalar**)array);CHKERRQ(ierr);
55     PetscFunctionReturn(0);
56   }
57   ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr);
58   ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gxm,&gym,&gzm);CHKERRQ(ierr);
59   ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr);
60 
61   /* Handle case where user passes in global vector as opposed to local */
62   ierr = VecGetLocalSize(vec,&N);CHKERRQ(ierr);
63   if (N == xm*ym*zm*dof) {
64     gxm = xm;
65     gym = ym;
66     gzm = zm;
67     gxs = xs;
68     gys = ys;
69     gzs = zs;
70   } else if (N != gxm*gym*gzm*dof) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vector local size %D is not compatible with DMDA local sizes %D %D\n",N,xm*ym*zm*dof,gxm*gym*gzm*dof);
71 
72   if (dim == 1) {
73     ierr = VecGetArray1d(vec,gxm*dof,gxs*dof,(PetscScalar**)array);CHKERRQ(ierr);
74   } else if (dim == 2) {
75     ierr = VecGetArray2d(vec,gym,gxm*dof,gys,gxs*dof,(PetscScalar***)array);CHKERRQ(ierr);
76   } else if (dim == 3) {
77     ierr = VecGetArray3d(vec,gzm,gym,gxm*dof,gzs,gys,gxs*dof,(PetscScalar****)array);CHKERRQ(ierr);
78   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim);
79   PetscFunctionReturn(0);
80 }
81 
82 #undef __FUNCT__
83 #define __FUNCT__ "DMDAVecRestoreArray"
84 /*@
85    DMDAVecRestoreArray - Restores a multiple dimension array obtained with DMDAVecGetArray()
86 
87    Logically collective on Vec
88 
89    Input Parameter:
90 +  da - the distributed array
91 .  vec - the vector, either a vector the same size as one obtained with
92          DMCreateGlobalVector() or DMCreateLocalVector()
93 -  array - the array, non-NULL pointer is zeroed
94 
95   Level: intermediate
96 
97   Fortran Notes: From Fortran use DMDAVecRestoreArrayF90()
98 
99 .keywords: distributed array, get, corners, nodes, local indices, coordinates
100 
101 .seealso: DMDAGetGhostCorners(), DMDAGetCorners(), VecGetArray(), VecRestoreArray(), DMDAVecGetArray()
102 @*/
103 PetscErrorCode  DMDAVecRestoreArray(DM da,Vec vec,void *array)
104 {
105   PetscErrorCode ierr;
106   PetscInt       xs,ys,zs,xm,ym,zm,gxs,gys,gzs,gxm,gym,gzm,N,dim,dof;
107 
108   PetscFunctionBegin;
109   PetscValidHeaderSpecific(da, DM_CLASSID, 1);
110   PetscValidHeaderSpecific(vec, VEC_CLASSID, 2);
111   PetscValidPointer(array, 3);
112   if (da->defaultSection) {
113     ierr = VecRestoreArray(vec,(PetscScalar**)array);CHKERRQ(ierr);
114     PetscFunctionReturn(0);
115   }
116   ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr);
117   ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gxm,&gym,&gzm);CHKERRQ(ierr);
118   ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr);
119 
120   /* Handle case where user passes in global vector as opposed to local */
121   ierr = VecGetLocalSize(vec,&N);CHKERRQ(ierr);
122   if (N == xm*ym*zm*dof) {
123     gxm = xm;
124     gym = ym;
125     gzm = zm;
126     gxs = xs;
127     gys = ys;
128     gzs = zs;
129   } else if (N != gxm*gym*gzm*dof) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vector local size %D is not compatible with DMDA local sizes %D %D\n",N,xm*ym*zm*dof,gxm*gym*gzm*dof);
130 
131   if (dim == 1) {
132     ierr = VecRestoreArray1d(vec,gxm*dof,gxs*dof,(PetscScalar**)array);CHKERRQ(ierr);
133   } else if (dim == 2) {
134     ierr = VecRestoreArray2d(vec,gym,gxm*dof,gys,gxs*dof,(PetscScalar***)array);CHKERRQ(ierr);
135   } else if (dim == 3) {
136     ierr = VecRestoreArray3d(vec,gzm,gym,gxm*dof,gzs,gys,gxs*dof,(PetscScalar****)array);CHKERRQ(ierr);
137   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim);
138   PetscFunctionReturn(0);
139 }
140 
141 #undef __FUNCT__
142 #define __FUNCT__ "DMDAVecGetArrayDOF"
143 /*@C
144    DMDAVecGetArrayDOF - Returns a multiple dimension array that shares data with
145       the underlying vector and is indexed using the global dimensions.
146 
147    Logically collective
148 
149    Input Parameter:
150 +  da - the distributed array
151 -  vec - the vector, either a vector the same size as one obtained with
152          DMCreateGlobalVector() or DMCreateLocalVector()
153 
154    Output Parameter:
155 .  array - the array
156 
157    Notes:
158     Call DMDAVecRestoreArrayDOF() once you have finished accessing the vector entries.
159 
160     In C, the indexing is "backwards" from what expects: array[k][j][i][DOF] NOT array[i][j][k][DOF]!
161 
162     In Fortran 90 you do not need a version of DMDAVecRestoreArrayDOF() just use  DMDAVecRestoreArrayF90() and declare your array with one higher dimension,
163     see src/dm/examples/tutorials/ex11f90.F
164 
165   Level: intermediate
166 
167 .keywords: distributed array, get, corners, nodes, local indices, coordinates
168 
169 .seealso: DMDAGetGhostCorners(), DMDAGetCorners(), VecGetArray(), VecRestoreArray(), DMDAVecRestoreArray(), DMDAVecGetArray(), DMDAVecRestoreArrayDOF()
170 @*/
171 PetscErrorCode  DMDAVecGetArrayDOF(DM da,Vec vec,void *array)
172 {
173   PetscErrorCode ierr;
174   PetscInt       xs,ys,zs,xm,ym,zm,gxs,gys,gzs,gxm,gym,gzm,N,dim,dof;
175 
176   PetscFunctionBegin;
177   ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr);
178   ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gxm,&gym,&gzm);CHKERRQ(ierr);
179   ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr);
180 
181   /* Handle case where user passes in global vector as opposed to local */
182   ierr = VecGetLocalSize(vec,&N);CHKERRQ(ierr);
183   if (N == xm*ym*zm*dof) {
184     gxm = xm;
185     gym = ym;
186     gzm = zm;
187     gxs = xs;
188     gys = ys;
189     gzs = zs;
190   } else if (N != gxm*gym*gzm*dof) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vector local size %D is not compatible with DMDA local sizes %D %D\n",N,xm*ym*zm*dof,gxm*gym*gzm*dof);
191 
192   if (dim == 1) {
193     ierr = VecGetArray2d(vec,gxm,dof,gxs,0,(PetscScalar***)array);CHKERRQ(ierr);
194   } else if (dim == 2) {
195     ierr = VecGetArray3d(vec,gym,gxm,dof,gys,gxs,0,(PetscScalar****)array);CHKERRQ(ierr);
196   } else if (dim == 3) {
197     ierr = VecGetArray4d(vec,gzm,gym,gxm,dof,gzs,gys,gxs,0,(PetscScalar*****)array);CHKERRQ(ierr);
198   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim);
199   PetscFunctionReturn(0);
200 }
201 
202 #undef __FUNCT__
203 #define __FUNCT__ "DMDAVecRestoreArrayDOF"
204 /*@
205    DMDAVecRestoreArrayDOF - Restores a multiple dimension array obtained with DMDAVecGetArrayDOF()
206 
207    Logically collective
208 
209    Input Parameter:
210 +  da - the distributed array
211 .  vec - the vector, either a vector the same size as one obtained with
212          DMCreateGlobalVector() or DMCreateLocalVector()
213 -  array - the array
214 
215   Level: intermediate
216 
217 .keywords: distributed array, get, corners, nodes, local indices, coordinates
218 
219 .seealso: DMDAGetGhostCorners(), DMDAGetCorners(), VecGetArray(), VecRestoreArray(), DMDAVecGetArray(), DMDAVecGetArrayDOF(), DMDAVecRestoreArrayDOF()
220 @*/
221 PetscErrorCode  DMDAVecRestoreArrayDOF(DM da,Vec vec,void *array)
222 {
223   PetscErrorCode ierr;
224   PetscInt       xs,ys,zs,xm,ym,zm,gxs,gys,gzs,gxm,gym,gzm,N,dim,dof;
225 
226   PetscFunctionBegin;
227   ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr);
228   ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gxm,&gym,&gzm);CHKERRQ(ierr);
229   ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr);
230 
231   /* Handle case where user passes in global vector as opposed to local */
232   ierr = VecGetLocalSize(vec,&N);CHKERRQ(ierr);
233   if (N == xm*ym*zm*dof) {
234     gxm = xm;
235     gym = ym;
236     gzm = zm;
237     gxs = xs;
238     gys = ys;
239     gzs = zs;
240   }
241 
242   if (dim == 1) {
243     ierr = VecRestoreArray2d(vec,gxm,dof,gxs,0,(PetscScalar***)array);CHKERRQ(ierr);
244   } else if (dim == 2) {
245     ierr = VecRestoreArray3d(vec,gym,gxm,dof,gys,gxs,0,(PetscScalar****)array);CHKERRQ(ierr);
246   } else if (dim == 3) {
247     ierr = VecRestoreArray4d(vec,gzm,gym,gxm,dof,gzs,gys,gxs,0,(PetscScalar*****)array);CHKERRQ(ierr);
248   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim);
249   PetscFunctionReturn(0);
250 }
251 
252 #undef __FUNCT__
253 #define __FUNCT__ "DMDAVecGetArrayRead"
254 /*@C
255    DMDAVecGetArrayRead - Returns a multiple dimension array that shares data with
256       the underlying vector and is indexed using the global dimensions.
257 
258    Not collective
259 
260    Input Parameter:
261 +  da - the distributed array
262 -  vec - the vector, either a vector the same size as one obtained with DMCreateGlobalVector() or DMCreateLocalVector()
263 
264    Output Parameter:
265 .  array - the array
266 
267    Notes:
268     Call DMDAVecRestoreArrayRead() once you have finished accessing the vector entries.
269 
270     In C, the indexing is "backwards" from what expects: array[k][j][i] NOT array[i][j][k]!
271 
272     If vec is a local vector (obtained with DMCreateLocalVector() etc) then the ghost point locations are accessible. If it is
273     a global vector then the ghost points are not accessible. Of course with the local vector you will have had to do the
274 
275     appropriate DMGlobalToLocalBegin() and DMGlobalToLocalEnd() to have correct values in the ghost locations.
276 
277   Fortran Notes: From Fortran use DMDAVecGetArrayReadF90() and pass for the array type PetscScalar,pointer :: array(:,...,:) of the appropriate
278        dimension. For a DMDA created with a dof of 1 use the dimension of the DMDA, for a DMDA created with a dof greater than 1 use one more than the
279        dimension of the DMDA. The order of the indices is array(xs:xs+xm-1,ys:ys+ym-1,zs:zs+zm-1) (when dof is 1) otherwise
280        array(0:dof-1,xs:xs+xm-1,ys:ys+ym-1,zs:zs+zm-1) where the values are obtained from
281        DMDAGetCorners() for a global array or DMDAGetGhostCorners() for a local array. Include petsc/finclude/petscdmda.h90 to access this routine.
282 
283   Due to bugs in the compiler DMDAVecGetArrayReadF90() does not work with gfortran versions before 4.5
284 
285   Level: intermediate
286 
287 .keywords: distributed array, get, corners, nodes, local indices, coordinates
288 
289 .seealso: DMDAGetGhostCorners(), DMDAGetCorners(), VecGetArray(), VecRestoreArray(), DMDAVecRestoreArray(), DMDAVecRestoreArrayDOF()
290           DMDAVecGetArrayDOF()
291 @*/
292 PetscErrorCode  DMDAVecGetArrayRead(DM da,Vec vec,void *array)
293 {
294   PetscErrorCode ierr;
295   PetscInt       xs,ys,zs,xm,ym,zm,gxs,gys,gzs,gxm,gym,gzm,N,dim,dof;
296 
297   PetscFunctionBegin;
298   PetscValidHeaderSpecific(da, DM_CLASSID, 1);
299   PetscValidHeaderSpecific(vec, VEC_CLASSID, 2);
300   PetscValidPointer(array, 3);
301   if (da->defaultSection) {
302     ierr = VecGetArrayRead(vec,(const PetscScalar**)array);CHKERRQ(ierr);
303     PetscFunctionReturn(0);
304   }
305   ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr);
306   ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gxm,&gym,&gzm);CHKERRQ(ierr);
307   ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr);
308 
309   /* Handle case where user passes in global vector as opposed to local */
310   ierr = VecGetLocalSize(vec,&N);CHKERRQ(ierr);
311   if (N == xm*ym*zm*dof) {
312     gxm = xm;
313     gym = ym;
314     gzm = zm;
315     gxs = xs;
316     gys = ys;
317     gzs = zs;
318   } else if (N != gxm*gym*gzm*dof) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vector local size %D is not compatible with DMDA local sizes %D %D\n",N,xm*ym*zm*dof,gxm*gym*gzm*dof);
319 
320   if (dim == 1) {
321     ierr = VecGetArray1dRead(vec,gxm*dof,gxs*dof,(PetscScalar**)array);CHKERRQ(ierr);
322   } else if (dim == 2) {
323     ierr = VecGetArray2dRead(vec,gym,gxm*dof,gys,gxs*dof,(PetscScalar***)array);CHKERRQ(ierr);
324   } else if (dim == 3) {
325     ierr = VecGetArray3dRead(vec,gzm,gym,gxm*dof,gzs,gys,gxs*dof,(PetscScalar****)array);CHKERRQ(ierr);
326   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim);
327   PetscFunctionReturn(0);
328 }
329 
330 #undef __FUNCT__
331 #define __FUNCT__ "DMDAVecRestoreArrayRead"
332 /*@
333    DMDAVecRestoreArrayRead - Restores a multiple dimension array obtained with DMDAVecGetArrayRead()
334 
335    Not collective
336 
337    Input Parameter:
338 +  da - the distributed array
339 .  vec - the vector, either a vector the same size as one obtained with
340          DMCreateGlobalVector() or DMCreateLocalVector()
341 -  array - the array, non-NULL pointer is zeroed
342 
343   Level: intermediate
344 
345   Fortran Notes: From Fortran use DMDAVecRestoreArrayReadF90()
346 
347 .keywords: distributed array, get, corners, nodes, local indices, coordinates
348 
349 .seealso: DMDAGetGhostCorners(), DMDAGetCorners(), VecGetArray(), VecRestoreArray(), DMDAVecGetArray()
350 @*/
351 PetscErrorCode  DMDAVecRestoreArrayRead(DM da,Vec vec,void *array)
352 {
353   PetscErrorCode ierr;
354   PetscInt       xs,ys,zs,xm,ym,zm,gxs,gys,gzs,gxm,gym,gzm,N,dim,dof;
355 
356   PetscFunctionBegin;
357   PetscValidHeaderSpecific(da, DM_CLASSID, 1);
358   PetscValidHeaderSpecific(vec, VEC_CLASSID, 2);
359   PetscValidPointer(array, 3);
360   if (da->defaultSection) {
361     ierr = VecRestoreArrayRead(vec,(const PetscScalar**)array);CHKERRQ(ierr);
362     PetscFunctionReturn(0);
363   }
364   ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr);
365   ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gxm,&gym,&gzm);CHKERRQ(ierr);
366   ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr);
367 
368   /* Handle case where user passes in global vector as opposed to local */
369   ierr = VecGetLocalSize(vec,&N);CHKERRQ(ierr);
370   if (N == xm*ym*zm*dof) {
371     gxm = xm;
372     gym = ym;
373     gzm = zm;
374     gxs = xs;
375     gys = ys;
376     gzs = zs;
377   } else if (N != gxm*gym*gzm*dof) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vector local size %D is not compatible with DMDA local sizes %D %D\n",N,xm*ym*zm*dof,gxm*gym*gzm*dof);
378 
379   if (dim == 1) {
380     ierr = VecRestoreArray1dRead(vec,gxm*dof,gxs*dof,(PetscScalar**)array);CHKERRQ(ierr);
381   } else if (dim == 2) {
382     ierr = VecRestoreArray2dRead(vec,gym,gxm*dof,gys,gxs*dof,(PetscScalar***)array);CHKERRQ(ierr);
383   } else if (dim == 3) {
384     ierr = VecRestoreArray3dRead(vec,gzm,gym,gxm*dof,gzs,gys,gxs*dof,(PetscScalar****)array);CHKERRQ(ierr);
385   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim);
386   PetscFunctionReturn(0);
387 }
388 
389 #undef __FUNCT__
390 #define __FUNCT__ "DMDAVecGetArrayDOFRead"
391 /*@C
392    DMDAVecGetArrayDOFRead - Returns a multiple dimension array that shares data with
393       the underlying vector and is indexed using the global dimensions.
394 
395    Not Collective
396 
397    Input Parameter:
398 +  da - the distributed array
399 -  vec - the vector, either a vector the same size as one obtained with
400          DMCreateGlobalVector() or DMCreateLocalVector()
401 
402    Output Parameter:
403 .  array - the array
404 
405    Notes:
406     Call DMDAVecRestoreArrayDOFRead() once you have finished accessing the vector entries.
407 
408     In C, the indexing is "backwards" from what expects: array[k][j][i][DOF] NOT array[i][j][k][DOF]!
409 
410     In Fortran 90 you do not need a version of DMDAVecRestoreArrayDOF() just use  DMDAVecRestoreArrayReadF90() and declare your array with one higher dimension,
411     see src/dm/examples/tutorials/ex11f90.F
412 
413   Level: intermediate
414 
415 .keywords: distributed array, get, corners, nodes, local indices, coordinates
416 
417 .seealso: DMDAGetGhostCorners(), DMDAGetCorners(), VecGetArray(), VecRestoreArray(), DMDAVecRestoreArray(), DMDAVecGetArray(), DMDAVecRestoreArrayDOF()
418 @*/
419 PetscErrorCode  DMDAVecGetArrayDOFRead(DM da,Vec vec,void *array)
420 {
421   PetscErrorCode ierr;
422   PetscInt       xs,ys,zs,xm,ym,zm,gxs,gys,gzs,gxm,gym,gzm,N,dim,dof;
423 
424   PetscFunctionBegin;
425   ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr);
426   ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gxm,&gym,&gzm);CHKERRQ(ierr);
427   ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr);
428 
429   /* Handle case where user passes in global vector as opposed to local */
430   ierr = VecGetLocalSize(vec,&N);CHKERRQ(ierr);
431   if (N == xm*ym*zm*dof) {
432     gxm = xm;
433     gym = ym;
434     gzm = zm;
435     gxs = xs;
436     gys = ys;
437     gzs = zs;
438   } else if (N != gxm*gym*gzm*dof) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vector local size %D is not compatible with DMDA local sizes %D %D\n",N,xm*ym*zm*dof,gxm*gym*gzm*dof);
439 
440   if (dim == 1) {
441     ierr = VecGetArray2dRead(vec,gxm,dof,gxs,0,(PetscScalar***)array);CHKERRQ(ierr);
442   } else if (dim == 2) {
443     ierr = VecGetArray3dRead(vec,gym,gxm,dof,gys,gxs,0,(PetscScalar****)array);CHKERRQ(ierr);
444   } else if (dim == 3) {
445     ierr = VecGetArray4dRead(vec,gzm,gym,gxm,dof,gzs,gys,gxs,0,(PetscScalar*****)array);CHKERRQ(ierr);
446   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim);
447   PetscFunctionReturn(0);
448 }
449 
450 #undef __FUNCT__
451 #define __FUNCT__ "DMDAVecRestoreArrayDOFRead"
452 /*@
453    DMDAVecRestoreArrayDOFRead - Restores a multiple dimension array obtained with DMDAVecGetArrayDOFRead()
454 
455    Not Collective
456 
457    Input Parameter:
458 +  da - the distributed array
459 .  vec - the vector, either a vector the same size as one obtained with
460          DMCreateGlobalVector() or DMCreateLocalVector()
461 -  array - the array
462 
463   Level: intermediate
464 
465 .keywords: distributed array, get, corners, nodes, local indices, coordinates
466 
467 .seealso: DMDAGetGhostCorners(), DMDAGetCorners(), VecGetArray(), VecRestoreArray(), DMDAVecGetArray(), DMDAVecGetArrayDOF(), DMDAVecRestoreArrayDOF()
468 @*/
469 PetscErrorCode  DMDAVecRestoreArrayDOFRead(DM da,Vec vec,void *array)
470 {
471   PetscErrorCode ierr;
472   PetscInt       xs,ys,zs,xm,ym,zm,gxs,gys,gzs,gxm,gym,gzm,N,dim,dof;
473 
474   PetscFunctionBegin;
475   ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr);
476   ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gxm,&gym,&gzm);CHKERRQ(ierr);
477   ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr);
478 
479   /* Handle case where user passes in global vector as opposed to local */
480   ierr = VecGetLocalSize(vec,&N);CHKERRQ(ierr);
481   if (N == xm*ym*zm*dof) {
482     gxm = xm;
483     gym = ym;
484     gzm = zm;
485     gxs = xs;
486     gys = ys;
487     gzs = zs;
488   }
489 
490   if (dim == 1) {
491     ierr = VecRestoreArray2dRead(vec,gxm,dof,gxs,0,(PetscScalar***)array);CHKERRQ(ierr);
492   } else if (dim == 2) {
493     ierr = VecRestoreArray3dRead(vec,gym,gxm,dof,gys,gxs,0,(PetscScalar****)array);CHKERRQ(ierr);
494   } else if (dim == 3) {
495     ierr = VecRestoreArray4dRead(vec,gzm,gym,gxm,dof,gzs,gys,gxs,0,(PetscScalar*****)array);CHKERRQ(ierr);
496   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim);
497   PetscFunctionReturn(0);
498 }
499 
500 
501 
502 
503 
504 
505 
506 
507 
508 
509 
510 
511 
512