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