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