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