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 */ 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 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 if (dd->nlocal < da->bind_below) {ierr = VecBindToCPU(*g,PETSC_TRUE);CHKERRQ(ierr);} 60 ierr = VecSetDM(*g, da);CHKERRQ(ierr); 61 #if defined(PETSC_HAVE_MATLAB_ENGINE) 62 if (dd->w == 1 && da->dim == 2) { 63 ierr = PetscObjectComposeFunction((PetscObject)*g,"PetscMatlabEnginePut_C",VecMatlabEnginePut_DA2d);CHKERRQ(ierr); 64 } 65 #endif 66 PetscFunctionReturn(0); 67 } 68 69 /*@ 70 DMDAGetNumCells - Get the number of cells in the local piece of the DMDA. This includes ghost cells. 71 72 Input Parameter: 73 . dm - The DM object 74 75 Output Parameters: 76 + numCellsX - The number of local cells in the x-direction 77 . numCellsY - The number of local cells in the y-direction 78 . numCellsZ - The number of local cells in the z-direction 79 - numCells - The number of local cells 80 81 Level: developer 82 83 .seealso: DMDAGetCellPoint() 84 @*/ 85 PetscErrorCode DMDAGetNumCells(DM dm, PetscInt *numCellsX, PetscInt *numCellsY, PetscInt *numCellsZ, PetscInt *numCells) 86 { 87 DM_DA *da = (DM_DA*) dm->data; 88 const PetscInt dim = dm->dim; 89 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 90 const PetscInt nC = (mx)*(dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1); 91 92 PetscFunctionBegin; 93 PetscValidHeaderSpecificType(dm, DM_CLASSID, 1,DMDA); 94 if (numCellsX) { 95 PetscValidIntPointer(numCellsX,2); 96 *numCellsX = mx; 97 } 98 if (numCellsY) { 99 PetscValidIntPointer(numCellsY,3); 100 *numCellsY = my; 101 } 102 if (numCellsZ) { 103 PetscValidIntPointer(numCellsZ,4); 104 *numCellsZ = mz; 105 } 106 if (numCells) { 107 PetscValidIntPointer(numCells,5); 108 *numCells = nC; 109 } 110 PetscFunctionReturn(0); 111 } 112 113 /*@ 114 DMDAGetCellPoint - Get the DM point corresponding to the tuple (i, j, k) in the DMDA 115 116 Input Parameters: 117 + dm - The DM object 118 - i,j,k - The global indices for the cell 119 120 Output Parameters: 121 . point - The local DM point 122 123 Level: developer 124 125 .seealso: DMDAGetNumCells() 126 @*/ 127 PetscErrorCode DMDAGetCellPoint(DM dm, PetscInt i, PetscInt j, PetscInt k, PetscInt *point) 128 { 129 const PetscInt dim = dm->dim; 130 DMDALocalInfo info; 131 PetscErrorCode ierr; 132 133 PetscFunctionBegin; 134 PetscValidHeaderSpecificType(dm, DM_CLASSID, 1,DMDA); 135 PetscValidIntPointer(point,5); 136 ierr = DMDAGetLocalInfo(dm, &info);CHKERRQ(ierr); 137 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);} 138 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);} 139 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);} 140 *point = i + (dim > 1 ? (j + (dim > 2 ? k*info.gym : 0))*info.gxm : 0); 141 PetscFunctionReturn(0); 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 PetscValidIntPointer(numVerticesX,2); 157 *numVerticesX = nVx; 158 } 159 if (numVerticesY) { 160 PetscValidIntPointer(numVerticesY,3); 161 *numVerticesY = nVy; 162 } 163 if (numVerticesZ) { 164 PetscValidIntPointer(numVerticesZ,4); 165 *numVerticesZ = nVz; 166 } 167 if (numVertices) { 168 PetscValidIntPointer(numVertices,5); 169 *numVertices = nV; 170 } 171 PetscFunctionReturn(0); 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 PetscValidIntPointer(numXFacesX,2); 189 *numXFacesX = nxF; 190 } 191 if (numXFaces) { 192 PetscValidIntPointer(numXFaces,3); 193 *numXFaces = nXF; 194 } 195 if (numYFacesY) { 196 PetscValidIntPointer(numYFacesY,4); 197 *numYFacesY = nyF; 198 } 199 if (numYFaces) { 200 PetscValidIntPointer(numYFaces,5); 201 *numYFaces = nYF; 202 } 203 if (numZFacesZ) { 204 PetscValidIntPointer(numZFacesZ,6); 205 *numZFacesZ = nzF; 206 } 207 if (numZFaces) { 208 PetscValidIntPointer(numZFaces,7); 209 *numZFaces = nZF; 210 } 211 PetscFunctionReturn(0); 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 PetscErrorCode ierr; 219 220 PetscFunctionBegin; 221 if (pStart) PetscValidIntPointer(pStart,3); 222 if (pEnd) PetscValidIntPointer(pEnd,4); 223 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); 224 ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr); 225 ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr); 226 if (height == 0) { 227 /* Cells */ 228 if (pStart) *pStart = 0; 229 if (pEnd) *pEnd = nC; 230 } else if (height == 1) { 231 /* Faces */ 232 if (pStart) *pStart = nC+nV; 233 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 234 } else if (height == dim) { 235 /* Vertices */ 236 if (pStart) *pStart = nC; 237 if (pEnd) *pEnd = nC+nV; 238 } else if (height < 0) { 239 /* All points */ 240 if (pStart) *pStart = 0; 241 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 242 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of height %d in the DA", height); 243 PetscFunctionReturn(0); 244 } 245 246 PetscErrorCode DMDAGetDepthStratum(DM dm, PetscInt depth, PetscInt *pStart, PetscInt *pEnd) 247 { 248 const PetscInt dim = dm->dim; 249 PetscInt nC, nV, nXF, nYF, nZF; 250 PetscErrorCode ierr; 251 252 PetscFunctionBegin; 253 if (pStart) PetscValidIntPointer(pStart,3); 254 if (pEnd) PetscValidIntPointer(pEnd,4); 255 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); 256 ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr); 257 ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr); 258 if (depth == dim) { 259 /* Cells */ 260 if (pStart) *pStart = 0; 261 if (pEnd) *pEnd = nC; 262 } else if (depth == dim-1) { 263 /* Faces */ 264 if (pStart) *pStart = nC+nV; 265 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 266 } else if (depth == 0) { 267 /* Vertices */ 268 if (pStart) *pStart = nC; 269 if (pEnd) *pEnd = nC+nV; 270 } else if (depth < 0) { 271 /* All points */ 272 if (pStart) *pStart = 0; 273 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 274 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of depth %d in the DA", depth); 275 PetscFunctionReturn(0); 276 } 277 278 PetscErrorCode DMDAGetConeSize(DM dm, PetscInt p, PetscInt *coneSize) 279 { 280 const PetscInt dim = dm->dim; 281 PetscInt nC, nV, nXF, nYF, nZF; 282 PetscErrorCode ierr; 283 284 PetscFunctionBegin; 285 *coneSize = 0; 286 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); 287 ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr); 288 ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr); 289 switch (dim) { 290 case 2: 291 if (p >= 0) { 292 if (p < nC) { 293 *coneSize = 4; 294 } else if (p < nC+nV) { 295 *coneSize = 0; 296 } else if (p < nC+nV+nXF+nYF+nZF) { 297 *coneSize = 2; 298 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF); 299 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p); 300 break; 301 case 3: 302 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D"); 303 } 304 PetscFunctionReturn(0); 305 } 306 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 } 349 PetscFunctionReturn(0); 350 } 351 352 PetscErrorCode DMDARestoreCone(DM dm, PetscInt p, PetscInt *cone[]) 353 { 354 PetscErrorCode ierr; 355 356 PetscFunctionBegin; 357 ierr = DMGetWorkArray(dm, 6, MPIU_INT, cone);CHKERRQ(ierr); 358 PetscFunctionReturn(0); 359 } 360 361 PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu) 362 { 363 DM_DA *da = (DM_DA *) dm->data; 364 Vec coordinates; 365 PetscSection section; 366 PetscScalar *coords; 367 PetscReal h[3]; 368 PetscInt dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k; 369 PetscErrorCode ierr; 370 371 PetscFunctionBegin; 372 PetscValidHeaderSpecificType(dm, DM_CLASSID, 1,DMDA); 373 ierr = DMDAGetInfo(dm, &dim, &M, &N, &P, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 374 if (dim > 3) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"The following code only works for dim <= 3"); 375 h[0] = (xu - xl)/M; 376 h[1] = (yu - yl)/N; 377 h[2] = (zu - zl)/P; 378 ierr = DMDAGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 379 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr); 380 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ion);CHKERRQ(ierr); 381 ierr = PetscSectionSetNumFields(section, 1);CHKERRQ(ierr); 382 ierr = PetscSectionSetFieldComponents(section, 0, dim);CHKERRQ(ierr); 383 ierr = PetscSectionSetChart(section, vStart, vEnd);CHKERRQ(ierr); 384 for (v = vStart; v < vEnd; ++v) { 385 ierr = PetscSectionSetDof(section, v, dim);CHKERRQ(ierr); 386 } 387 ierr = PetscSectionSetUp(section);CHKERRQ(ierr); 388 ierr = PetscSectionGetStorageSize(section, &size);CHKERRQ(ierr); 389 ierr = VecCreateSeq(PETSC_COMM_SELF, size, &coordinates);CHKERRQ(ierr); 390 ierr = PetscObjectSetName((PetscObject)coordinates,"coordinates");CHKERRQ(ierr); 391 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 392 for (k = 0; k < nVz; ++k) { 393 PetscInt ind[3], d, off; 394 395 ind[0] = 0; 396 ind[1] = 0; 397 ind[2] = k + da->zs; 398 for (j = 0; j < nVy; ++j) { 399 ind[1] = j + da->ys; 400 for (i = 0; i < nVx; ++i) { 401 const PetscInt vertex = (k*nVy + j)*nVx + i + vStart; 402 403 ierr = PetscSectionGetOffset(section, vertex, &off);CHKERRQ(ierr); 404 ind[0] = i + da->xs; 405 for (d = 0; d < dim; ++d) { 406 coords[off+d] = h[d]*ind[d]; 407 } 408 } 409 } 410 } 411 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 412 ierr = DMSetCoordinateSection(dm, PETSC_DETERMINE, section);CHKERRQ(ierr); 413 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 414 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr); 415 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 416 PetscFunctionReturn(0); 417 } 418 419 /* ------------------------------------------------------------------- */ 420 421 /*@C 422 DMDAGetArray - Gets a work array for a DMDA 423 424 Input Parameters: 425 + da - information about my local patch 426 - ghosted - do you want arrays for the ghosted or nonghosted patch 427 428 Output Parameters: 429 . vptr - array data structured 430 431 Note: The vector values are NOT initialized and may have garbage in them, so you may need 432 to zero them. 433 434 Level: advanced 435 436 .seealso: DMDARestoreArray() 437 438 @*/ 439 PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr) 440 { 441 PetscErrorCode ierr; 442 PetscInt j,i,xs,ys,xm,ym,zs,zm; 443 char *iarray_start; 444 void **iptr = (void**)vptr; 445 DM_DA *dd = (DM_DA*)da->data; 446 447 PetscFunctionBegin; 448 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA); 449 if (ghosted) { 450 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 451 if (dd->arrayghostedin[i]) { 452 *iptr = dd->arrayghostedin[i]; 453 iarray_start = (char*)dd->startghostedin[i]; 454 dd->arrayghostedin[i] = NULL; 455 dd->startghostedin[i] = NULL; 456 457 goto done; 458 } 459 } 460 xs = dd->Xs; 461 ys = dd->Ys; 462 zs = dd->Zs; 463 xm = dd->Xe-dd->Xs; 464 ym = dd->Ye-dd->Ys; 465 zm = dd->Ze-dd->Zs; 466 } else { 467 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 468 if (dd->arrayin[i]) { 469 *iptr = dd->arrayin[i]; 470 iarray_start = (char*)dd->startin[i]; 471 dd->arrayin[i] = NULL; 472 dd->startin[i] = NULL; 473 474 goto done; 475 } 476 } 477 xs = dd->xs; 478 ys = dd->ys; 479 zs = dd->zs; 480 xm = dd->xe-dd->xs; 481 ym = dd->ye-dd->ys; 482 zm = dd->ze-dd->zs; 483 } 484 485 switch (da->dim) { 486 case 1: { 487 void *ptr; 488 489 ierr = PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 490 491 ptr = (void*)(iarray_start - xs*sizeof(PetscScalar)); 492 *iptr = (void*)ptr; 493 break; 494 } 495 case 2: { 496 void **ptr; 497 498 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 499 500 ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*)); 501 for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs); 502 *iptr = (void*)ptr; 503 break; 504 } 505 case 3: { 506 void ***ptr,**bptr; 507 508 ierr = PetscMalloc((zm+1)*sizeof(void**)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 509 510 ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*)); 511 bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**)); 512 for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys); 513 for (i=zs; i<zs+zm; i++) { 514 for (j=ys; j<ys+ym; j++) { 515 ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs); 516 } 517 } 518 *iptr = (void*)ptr; 519 break; 520 } 521 default: 522 SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",da->dim); 523 } 524 525 done: 526 /* add arrays to the checked out list */ 527 if (ghosted) { 528 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 529 if (!dd->arrayghostedout[i]) { 530 dd->arrayghostedout[i] = *iptr; 531 dd->startghostedout[i] = iarray_start; 532 break; 533 } 534 } 535 } else { 536 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 537 if (!dd->arrayout[i]) { 538 dd->arrayout[i] = *iptr; 539 dd->startout[i] = iarray_start; 540 break; 541 } 542 } 543 } 544 PetscFunctionReturn(0); 545 } 546 547 /*@C 548 DMDARestoreArray - Restores an array of derivative types for a DMDA 549 550 Input Parameters: 551 + da - information about my local patch 552 . ghosted - do you want arrays for the ghosted or nonghosted patch 553 - vptr - array data structured 554 555 Level: advanced 556 557 .seealso: DMDAGetArray() 558 559 @*/ 560 PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr) 561 { 562 PetscInt i; 563 void **iptr = (void**)vptr,*iarray_start = NULL; 564 DM_DA *dd = (DM_DA*)da->data; 565 566 PetscFunctionBegin; 567 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA); 568 if (ghosted) { 569 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 570 if (dd->arrayghostedout[i] == *iptr) { 571 iarray_start = dd->startghostedout[i]; 572 dd->arrayghostedout[i] = NULL; 573 dd->startghostedout[i] = NULL; 574 break; 575 } 576 } 577 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 578 if (!dd->arrayghostedin[i]) { 579 dd->arrayghostedin[i] = *iptr; 580 dd->startghostedin[i] = iarray_start; 581 break; 582 } 583 } 584 } else { 585 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 586 if (dd->arrayout[i] == *iptr) { 587 iarray_start = dd->startout[i]; 588 dd->arrayout[i] = NULL; 589 dd->startout[i] = NULL; 590 break; 591 } 592 } 593 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 594 if (!dd->arrayin[i]) { 595 dd->arrayin[i] = *iptr; 596 dd->startin[i] = iarray_start; 597 break; 598 } 599 } 600 } 601 PetscFunctionReturn(0); 602 } 603 604