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