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