1
2 /*
3 Code for manipulating distributed regular arrays in parallel.
4 */
5
6 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/
7 #include <petscbt.h>
8 #include <petscsf.h>
9 #include <petscds.h>
10 #include <petscfe.h>
11
12 /*
13 This allows the DMDA vectors to properly tell MATLAB their dimensions
14 */
15 #if defined(PETSC_HAVE_MATLAB_ENGINE)
16 #include <engine.h> /* MATLAB include file */
17 #include <mex.h> /* MATLAB include file */
VecMatlabEnginePut_DA2d(PetscObject obj,void * mengine)18 static PetscErrorCode VecMatlabEnginePut_DA2d(PetscObject obj,void *mengine)
19 {
20 PetscErrorCode ierr;
21 PetscInt n,m;
22 Vec vec = (Vec)obj;
23 PetscScalar *array;
24 mxArray *mat;
25 DM da;
26
27 PetscFunctionBegin;
28 ierr = VecGetDM(vec, &da);CHKERRQ(ierr);
29 if (!da) SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_ARG_WRONGSTATE,"Vector not associated with a DMDA");
30 ierr = DMDAGetGhostCorners(da,0,0,0,&m,&n,0);CHKERRQ(ierr);
31
32 ierr = VecGetArray(vec,&array);CHKERRQ(ierr);
33 #if !defined(PETSC_USE_COMPLEX)
34 mat = mxCreateDoubleMatrix(m,n,mxREAL);
35 #else
36 mat = mxCreateDoubleMatrix(m,n,mxCOMPLEX);
37 #endif
38 ierr = PetscArraycpy(mxGetPr(mat),array,n*m);CHKERRQ(ierr);
39 ierr = PetscObjectName(obj);CHKERRQ(ierr);
40 engPutVariable((Engine*)mengine,obj->name,mat);
41
42 ierr = VecRestoreArray(vec,&array);CHKERRQ(ierr);
43 PetscFunctionReturn(0);
44 }
45 #endif
46
DMCreateLocalVector_DA(DM da,Vec * g)47 PetscErrorCode DMCreateLocalVector_DA(DM da,Vec *g)
48 {
49 PetscErrorCode ierr;
50 DM_DA *dd = (DM_DA*)da->data;
51
52 PetscFunctionBegin;
53 PetscValidHeaderSpecific(da,DM_CLASSID,1);
54 PetscValidPointer(g,2);
55 ierr = VecCreate(PETSC_COMM_SELF,g);CHKERRQ(ierr);
56 ierr = VecSetSizes(*g,dd->nlocal,PETSC_DETERMINE);CHKERRQ(ierr);
57 ierr = VecSetBlockSize(*g,dd->w);CHKERRQ(ierr);
58 ierr = VecSetType(*g,da->vectype);CHKERRQ(ierr);
59 ierr = VecSetDM(*g, da);CHKERRQ(ierr);
60 #if defined(PETSC_HAVE_MATLAB_ENGINE)
61 if (dd->w == 1 && da->dim == 2) {
62 ierr = PetscObjectComposeFunction((PetscObject)*g,"PetscMatlabEnginePut_C",VecMatlabEnginePut_DA2d);CHKERRQ(ierr);
63 }
64 #endif
65 PetscFunctionReturn(0);
66 }
67
68 /*@
69 DMDAGetNumCells - Get the number of cells in the local piece of the DMDA. This includes ghost cells.
70
71 Input Parameter:
72 . dm - The DM object
73
74 Output Parameters:
75 + numCellsX - The number of local cells in the x-direction
76 . numCellsY - The number of local cells in the y-direction
77 . numCellsZ - The number of local cells in the z-direction
78 - numCells - The number of local cells
79
80 Level: developer
81
82 .seealso: DMDAGetCellPoint()
83 @*/
DMDAGetNumCells(DM dm,PetscInt * numCellsX,PetscInt * numCellsY,PetscInt * numCellsZ,PetscInt * numCells)84 PetscErrorCode DMDAGetNumCells(DM dm, PetscInt *numCellsX, PetscInt *numCellsY, PetscInt *numCellsZ, PetscInt *numCells)
85 {
86 DM_DA *da = (DM_DA*) dm->data;
87 const PetscInt dim = dm->dim;
88 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
89 const PetscInt nC = (mx)*(dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1);
90
91 PetscFunctionBegin;
92 PetscValidHeaderSpecificType(dm, DM_CLASSID, 1,DMDA);
93 if (numCellsX) {
94 PetscValidIntPointer(numCellsX,2);
95 *numCellsX = mx;
96 }
97 if (numCellsY) {
98 PetscValidIntPointer(numCellsY,3);
99 *numCellsY = my;
100 }
101 if (numCellsZ) {
102 PetscValidIntPointer(numCellsZ,4);
103 *numCellsZ = mz;
104 }
105 if (numCells) {
106 PetscValidIntPointer(numCells,5);
107 *numCells = nC;
108 }
109 PetscFunctionReturn(0);
110 }
111
112 /*@
113 DMDAGetCellPoint - Get the DM point corresponding to the tuple (i, j, k) in the DMDA
114
115 Input Parameters:
116 + dm - The DM object
117 - i,j,k - The global indices for the cell
118
119 Output Parameters:
120 . point - The local DM point
121
122 Level: developer
123
124 .seealso: DMDAGetNumCells()
125 @*/
DMDAGetCellPoint(DM dm,PetscInt i,PetscInt j,PetscInt k,PetscInt * point)126 PetscErrorCode DMDAGetCellPoint(DM dm, PetscInt i, PetscInt j, PetscInt k, PetscInt *point)
127 {
128 const PetscInt dim = dm->dim;
129 DMDALocalInfo info;
130 PetscErrorCode ierr;
131
132 PetscFunctionBegin;
133 PetscValidHeaderSpecificType(dm, DM_CLASSID, 1,DMDA);
134 PetscValidIntPointer(point,5);
135 ierr = DMDAGetLocalInfo(dm, &info);CHKERRQ(ierr);
136 if (dim > 0) {if ((i < info.gxs) || (i >= info.gxs+info.gxm)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "X index %d not in [%d, %d)", i, info.gxs, info.gxs+info.gxm);}
137 if (dim > 1) {if ((j < info.gys) || (j >= info.gys+info.gym)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Y index %d not in [%d, %d)", j, info.gys, info.gys+info.gym);}
138 if (dim > 2) {if ((k < info.gzs) || (k >= info.gzs+info.gzm)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Z index %d not in [%d, %d)", k, info.gzs, info.gzs+info.gzm);}
139 *point = i + (dim > 1 ? (j + (dim > 2 ? k*info.gym : 0))*info.gxm : 0);
140 PetscFunctionReturn(0);
141 }
142
DMDAGetNumVertices(DM dm,PetscInt * numVerticesX,PetscInt * numVerticesY,PetscInt * numVerticesZ,PetscInt * numVertices)143 PetscErrorCode DMDAGetNumVertices(DM dm, PetscInt *numVerticesX, PetscInt *numVerticesY, PetscInt *numVerticesZ, PetscInt *numVertices)
144 {
145 DM_DA *da = (DM_DA*) dm->data;
146 const PetscInt dim = dm->dim;
147 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
148 const PetscInt nVx = mx+1;
149 const PetscInt nVy = dim > 1 ? (my+1) : 1;
150 const PetscInt nVz = dim > 2 ? (mz+1) : 1;
151 const PetscInt nV = nVx*nVy*nVz;
152
153 PetscFunctionBegin;
154 if (numVerticesX) {
155 PetscValidIntPointer(numVerticesX,2);
156 *numVerticesX = nVx;
157 }
158 if (numVerticesY) {
159 PetscValidIntPointer(numVerticesY,3);
160 *numVerticesY = nVy;
161 }
162 if (numVerticesZ) {
163 PetscValidIntPointer(numVerticesZ,4);
164 *numVerticesZ = nVz;
165 }
166 if (numVertices) {
167 PetscValidIntPointer(numVertices,5);
168 *numVertices = nV;
169 }
170 PetscFunctionReturn(0);
171 }
172
DMDAGetNumFaces(DM dm,PetscInt * numXFacesX,PetscInt * numXFaces,PetscInt * numYFacesY,PetscInt * numYFaces,PetscInt * numZFacesZ,PetscInt * numZFaces)173 PetscErrorCode DMDAGetNumFaces(DM dm, PetscInt *numXFacesX, PetscInt *numXFaces, PetscInt *numYFacesY, PetscInt *numYFaces, PetscInt *numZFacesZ, PetscInt *numZFaces)
174 {
175 DM_DA *da = (DM_DA*) dm->data;
176 const PetscInt dim = dm->dim;
177 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
178 const PetscInt nxF = (dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1);
179 const PetscInt nXF = (mx+1)*nxF;
180 const PetscInt nyF = mx*(dim > 2 ? mz : 1);
181 const PetscInt nYF = dim > 1 ? (my+1)*nyF : 0;
182 const PetscInt nzF = mx*(dim > 1 ? my : 0);
183 const PetscInt nZF = dim > 2 ? (mz+1)*nzF : 0;
184
185 PetscFunctionBegin;
186 if (numXFacesX) {
187 PetscValidIntPointer(numXFacesX,2);
188 *numXFacesX = nxF;
189 }
190 if (numXFaces) {
191 PetscValidIntPointer(numXFaces,3);
192 *numXFaces = nXF;
193 }
194 if (numYFacesY) {
195 PetscValidIntPointer(numYFacesY,4);
196 *numYFacesY = nyF;
197 }
198 if (numYFaces) {
199 PetscValidIntPointer(numYFaces,5);
200 *numYFaces = nYF;
201 }
202 if (numZFacesZ) {
203 PetscValidIntPointer(numZFacesZ,6);
204 *numZFacesZ = nzF;
205 }
206 if (numZFaces) {
207 PetscValidIntPointer(numZFaces,7);
208 *numZFaces = nZF;
209 }
210 PetscFunctionReturn(0);
211 }
212
DMDAGetHeightStratum(DM dm,PetscInt height,PetscInt * pStart,PetscInt * pEnd)213 PetscErrorCode DMDAGetHeightStratum(DM dm, PetscInt height, PetscInt *pStart, PetscInt *pEnd)
214 {
215 const PetscInt dim = dm->dim;
216 PetscInt nC, nV, nXF, nYF, nZF;
217 PetscErrorCode ierr;
218
219 PetscFunctionBegin;
220 if (pStart) PetscValidIntPointer(pStart,3);
221 if (pEnd) PetscValidIntPointer(pEnd,4);
222 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr);
223 ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr);
224 ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr);
225 if (height == 0) {
226 /* Cells */
227 if (pStart) *pStart = 0;
228 if (pEnd) *pEnd = nC;
229 } else if (height == 1) {
230 /* Faces */
231 if (pStart) *pStart = nC+nV;
232 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
233 } else if (height == dim) {
234 /* Vertices */
235 if (pStart) *pStart = nC;
236 if (pEnd) *pEnd = nC+nV;
237 } else if (height < 0) {
238 /* All points */
239 if (pStart) *pStart = 0;
240 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
241 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of height %d in the DA", height);
242 PetscFunctionReturn(0);
243 }
244
DMDAGetDepthStratum(DM dm,PetscInt depth,PetscInt * pStart,PetscInt * pEnd)245 PetscErrorCode DMDAGetDepthStratum(DM dm, PetscInt depth, PetscInt *pStart, PetscInt *pEnd)
246 {
247 const PetscInt dim = dm->dim;
248 PetscInt nC, nV, nXF, nYF, nZF;
249 PetscErrorCode ierr;
250
251 PetscFunctionBegin;
252 if (pStart) PetscValidIntPointer(pStart,3);
253 if (pEnd) PetscValidIntPointer(pEnd,4);
254 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr);
255 ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr);
256 ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr);
257 if (depth == dim) {
258 /* Cells */
259 if (pStart) *pStart = 0;
260 if (pEnd) *pEnd = nC;
261 } else if (depth == dim-1) {
262 /* Faces */
263 if (pStart) *pStart = nC+nV;
264 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
265 } else if (depth == 0) {
266 /* Vertices */
267 if (pStart) *pStart = nC;
268 if (pEnd) *pEnd = nC+nV;
269 } else if (depth < 0) {
270 /* All points */
271 if (pStart) *pStart = 0;
272 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
273 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of depth %d in the DA", depth);
274 PetscFunctionReturn(0);
275 }
276
DMDAGetConeSize(DM dm,PetscInt p,PetscInt * coneSize)277 PetscErrorCode DMDAGetConeSize(DM dm, PetscInt p, PetscInt *coneSize)
278 {
279 const PetscInt dim = dm->dim;
280 PetscInt nC, nV, nXF, nYF, nZF;
281 PetscErrorCode ierr;
282
283 PetscFunctionBegin;
284 *coneSize = 0;
285 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr);
286 ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr);
287 ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr);
288 switch (dim) {
289 case 2:
290 if (p >= 0) {
291 if (p < nC) {
292 *coneSize = 4;
293 } else if (p < nC+nV) {
294 *coneSize = 0;
295 } else if (p < nC+nV+nXF+nYF+nZF) {
296 *coneSize = 2;
297 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF);
298 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p);
299 break;
300 case 3:
301 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
302 break;
303 }
304 PetscFunctionReturn(0);
305 }
306
DMDAGetCone(DM dm,PetscInt p,PetscInt * cone[])307 PetscErrorCode DMDAGetCone(DM dm, PetscInt p, PetscInt *cone[])
308 {
309 const PetscInt dim = dm->dim;
310 PetscInt nCx, nCy, nCz, nC, nVx, nVy, nVz, nV, nxF, nyF, nzF, nXF, nYF, nZF;
311 PetscErrorCode ierr;
312
313 PetscFunctionBegin;
314 if (!cone) {ierr = DMGetWorkArray(dm, 6, MPIU_INT, cone);CHKERRQ(ierr);}
315 ierr = DMDAGetNumCells(dm, &nCx, &nCy, &nCz, &nC);CHKERRQ(ierr);
316 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr);
317 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr);
318 switch (dim) {
319 case 2:
320 if (p >= 0) {
321 if (p < nC) {
322 const PetscInt cy = p / nCx;
323 const PetscInt cx = p % nCx;
324
325 (*cone)[0] = cy*nxF + cx + nC+nV;
326 (*cone)[1] = cx*nyF + cy + nyF + nC+nV+nXF;
327 (*cone)[2] = cy*nxF + cx + nxF + nC+nV;
328 (*cone)[3] = cx*nyF + cy + nC+nV+nXF;
329 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do cell cones");
330 } else if (p < nC+nV) {
331 } else if (p < nC+nV+nXF) {
332 const PetscInt fy = (p - nC+nV) / nxF;
333 const PetscInt fx = (p - nC+nV) % nxF;
334
335 (*cone)[0] = fy*nVx + fx + nC;
336 (*cone)[1] = fy*nVx + fx + 1 + nC;
337 } else if (p < nC+nV+nXF+nYF) {
338 const PetscInt fx = (p - nC+nV+nXF) / nyF;
339 const PetscInt fy = (p - nC+nV+nXF) % nyF;
340
341 (*cone)[0] = fy*nVx + fx + nC;
342 (*cone)[1] = fy*nVx + fx + nVx + nC;
343 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF);
344 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p);
345 break;
346 case 3:
347 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
348 break;
349 }
350 PetscFunctionReturn(0);
351 }
352
DMDARestoreCone(DM dm,PetscInt p,PetscInt * cone[])353 PetscErrorCode DMDARestoreCone(DM dm, PetscInt p, PetscInt *cone[])
354 {
355 PetscErrorCode ierr;
356
357 PetscFunctionBegin;
358 ierr = DMGetWorkArray(dm, 6, MPIU_INT, cone);CHKERRQ(ierr);
359 PetscFunctionReturn(0);
360 }
361
DMDASetVertexCoordinates(DM dm,PetscReal xl,PetscReal xu,PetscReal yl,PetscReal yu,PetscReal zl,PetscReal zu)362 PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu)
363 {
364 DM_DA *da = (DM_DA *) dm->data;
365 Vec coordinates;
366 PetscSection section;
367 PetscScalar *coords;
368 PetscReal h[3];
369 PetscInt dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k;
370 PetscErrorCode ierr;
371
372 PetscFunctionBegin;
373 PetscValidHeaderSpecificType(dm, DM_CLASSID, 1,DMDA);
374 ierr = DMDAGetInfo(dm, &dim, &M, &N, &P, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
375 if (dim > 3) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"The following code only works for dim <= 3");
376 h[0] = (xu - xl)/M;
377 h[1] = (yu - yl)/N;
378 h[2] = (zu - zl)/P;
379 ierr = DMDAGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
380 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr);
381 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ion);CHKERRQ(ierr);
382 ierr = PetscSectionSetNumFields(section, 1);CHKERRQ(ierr);
383 ierr = PetscSectionSetFieldComponents(section, 0, dim);CHKERRQ(ierr);
384 ierr = PetscSectionSetChart(section, vStart, vEnd);CHKERRQ(ierr);
385 for (v = vStart; v < vEnd; ++v) {
386 ierr = PetscSectionSetDof(section, v, dim);CHKERRQ(ierr);
387 }
388 ierr = PetscSectionSetUp(section);CHKERRQ(ierr);
389 ierr = PetscSectionGetStorageSize(section, &size);CHKERRQ(ierr);
390 ierr = VecCreateSeq(PETSC_COMM_SELF, size, &coordinates);CHKERRQ(ierr);
391 ierr = PetscObjectSetName((PetscObject)coordinates,"coordinates");CHKERRQ(ierr);
392 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
393 for (k = 0; k < nVz; ++k) {
394 PetscInt ind[3], d, off;
395
396 ind[0] = 0;
397 ind[1] = 0;
398 ind[2] = k + da->zs;
399 for (j = 0; j < nVy; ++j) {
400 ind[1] = j + da->ys;
401 for (i = 0; i < nVx; ++i) {
402 const PetscInt vertex = (k*nVy + j)*nVx + i + vStart;
403
404 ierr = PetscSectionGetOffset(section, vertex, &off);CHKERRQ(ierr);
405 ind[0] = i + da->xs;
406 for (d = 0; d < dim; ++d) {
407 coords[off+d] = h[d]*ind[d];
408 }
409 }
410 }
411 }
412 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
413 ierr = DMSetCoordinateSection(dm, PETSC_DETERMINE, section);CHKERRQ(ierr);
414 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr);
415 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr);
416 ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
417 PetscFunctionReturn(0);
418 }
419
420 /* ------------------------------------------------------------------- */
421
422 /*@C
423 DMDAGetArray - Gets a work array for a DMDA
424
425 Input Parameter:
426 + da - information about my local patch
427 - ghosted - do you want arrays for the ghosted or nonghosted patch
428
429 Output Parameters:
430 . vptr - array data structured
431
432 Note: The vector values are NOT initialized and may have garbage in them, so you may need
433 to zero them.
434
435 Level: advanced
436
437 .seealso: DMDARestoreArray()
438
439 @*/
DMDAGetArray(DM da,PetscBool ghosted,void * vptr)440 PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr)
441 {
442 PetscErrorCode ierr;
443 PetscInt j,i,xs,ys,xm,ym,zs,zm;
444 char *iarray_start;
445 void **iptr = (void**)vptr;
446 DM_DA *dd = (DM_DA*)da->data;
447
448 PetscFunctionBegin;
449 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
450 if (ghosted) {
451 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
452 if (dd->arrayghostedin[i]) {
453 *iptr = dd->arrayghostedin[i];
454 iarray_start = (char*)dd->startghostedin[i];
455 dd->arrayghostedin[i] = NULL;
456 dd->startghostedin[i] = NULL;
457
458 goto done;
459 }
460 }
461 xs = dd->Xs;
462 ys = dd->Ys;
463 zs = dd->Zs;
464 xm = dd->Xe-dd->Xs;
465 ym = dd->Ye-dd->Ys;
466 zm = dd->Ze-dd->Zs;
467 } else {
468 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
469 if (dd->arrayin[i]) {
470 *iptr = dd->arrayin[i];
471 iarray_start = (char*)dd->startin[i];
472 dd->arrayin[i] = NULL;
473 dd->startin[i] = NULL;
474
475 goto done;
476 }
477 }
478 xs = dd->xs;
479 ys = dd->ys;
480 zs = dd->zs;
481 xm = dd->xe-dd->xs;
482 ym = dd->ye-dd->ys;
483 zm = dd->ze-dd->zs;
484 }
485
486 switch (da->dim) {
487 case 1: {
488 void *ptr;
489
490 ierr = PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr);
491
492 ptr = (void*)(iarray_start - xs*sizeof(PetscScalar));
493 *iptr = (void*)ptr;
494 break;
495 }
496 case 2: {
497 void **ptr;
498
499 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr);
500
501 ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*));
502 for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs);
503 *iptr = (void*)ptr;
504 break;
505 }
506 case 3: {
507 void ***ptr,**bptr;
508
509 ierr = PetscMalloc((zm+1)*sizeof(void**)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr);
510
511 ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*));
512 bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**));
513 for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys);
514 for (i=zs; i<zs+zm; i++) {
515 for (j=ys; j<ys+ym; j++) {
516 ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs);
517 }
518 }
519 *iptr = (void*)ptr;
520 break;
521 }
522 default:
523 SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",da->dim);
524 }
525
526 done:
527 /* add arrays to the checked out list */
528 if (ghosted) {
529 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
530 if (!dd->arrayghostedout[i]) {
531 dd->arrayghostedout[i] = *iptr;
532 dd->startghostedout[i] = iarray_start;
533 break;
534 }
535 }
536 } else {
537 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
538 if (!dd->arrayout[i]) {
539 dd->arrayout[i] = *iptr;
540 dd->startout[i] = iarray_start;
541 break;
542 }
543 }
544 }
545 PetscFunctionReturn(0);
546 }
547
548 /*@C
549 DMDARestoreArray - Restores an array of derivative types for a DMDA
550
551 Input Parameter:
552 + da - information about my local patch
553 . ghosted - do you want arrays for the ghosted or nonghosted patch
554 - vptr - array data structured
555
556 Level: advanced
557
558 .seealso: DMDAGetArray()
559
560 @*/
DMDARestoreArray(DM da,PetscBool ghosted,void * vptr)561 PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr)
562 {
563 PetscInt i;
564 void **iptr = (void**)vptr,*iarray_start = NULL;
565 DM_DA *dd = (DM_DA*)da->data;
566
567 PetscFunctionBegin;
568 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
569 if (ghosted) {
570 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
571 if (dd->arrayghostedout[i] == *iptr) {
572 iarray_start = dd->startghostedout[i];
573 dd->arrayghostedout[i] = NULL;
574 dd->startghostedout[i] = NULL;
575 break;
576 }
577 }
578 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
579 if (!dd->arrayghostedin[i]) {
580 dd->arrayghostedin[i] = *iptr;
581 dd->startghostedin[i] = iarray_start;
582 break;
583 }
584 }
585 } else {
586 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
587 if (dd->arrayout[i] == *iptr) {
588 iarray_start = dd->startout[i];
589 dd->arrayout[i] = NULL;
590 dd->startout[i] = NULL;
591 break;
592 }
593 }
594 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
595 if (!dd->arrayin[i]) {
596 dd->arrayin[i] = *iptr;
597 dd->startin[i] = iarray_start;
598 break;
599 }
600 }
601 }
602 PetscFunctionReturn(0);
603 }
604
605