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 <petscsf.h> 8 9 /* 10 This allows the DMDA vectors to properly tell MATLAB their dimensions 11 */ 12 #if defined(PETSC_HAVE_MATLAB_ENGINE) 13 #include <engine.h> /* MATLAB include file */ 14 #include <mex.h> /* MATLAB include file */ 15 #undef __FUNCT__ 16 #define __FUNCT__ "VecMatlabEnginePut_DA2d" 17 static PetscErrorCode VecMatlabEnginePut_DA2d(PetscObject obj,void *mengine) 18 { 19 PetscErrorCode ierr; 20 PetscInt n,m; 21 Vec vec = (Vec)obj; 22 PetscScalar *array; 23 mxArray *mat; 24 DM da; 25 26 PetscFunctionBegin; 27 ierr = VecGetDM(vec, &da);CHKERRQ(ierr); 28 if (!da) SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_ARG_WRONGSTATE,"Vector not associated with a DMDA"); 29 ierr = DMDAGetGhostCorners(da,0,0,0,&m,&n,0);CHKERRQ(ierr); 30 31 ierr = VecGetArray(vec,&array);CHKERRQ(ierr); 32 #if !defined(PETSC_USE_COMPLEX) 33 mat = mxCreateDoubleMatrix(m,n,mxREAL); 34 #else 35 mat = mxCreateDoubleMatrix(m,n,mxCOMPLEX); 36 #endif 37 ierr = PetscMemcpy(mxGetPr(mat),array,n*m*sizeof(PetscScalar));CHKERRQ(ierr); 38 ierr = PetscObjectName(obj);CHKERRQ(ierr); 39 engPutVariable((Engine*)mengine,obj->name,mat); 40 41 ierr = VecRestoreArray(vec,&array);CHKERRQ(ierr); 42 PetscFunctionReturn(0); 43 } 44 #endif 45 46 47 #undef __FUNCT__ 48 #define __FUNCT__ "DMCreateLocalVector_DA" 49 PetscErrorCode DMCreateLocalVector_DA(DM da,Vec *g) 50 { 51 PetscErrorCode ierr; 52 DM_DA *dd = (DM_DA*)da->data; 53 54 PetscFunctionBegin; 55 PetscValidHeaderSpecific(da,DM_CLASSID,1); 56 PetscValidPointer(g,2); 57 if (da->defaultSection) { 58 ierr = DMCreateLocalVector_Section_Private(da,g);CHKERRQ(ierr); 59 } else { 60 ierr = VecCreate(PETSC_COMM_SELF,g);CHKERRQ(ierr); 61 ierr = VecSetSizes(*g,dd->nlocal,PETSC_DETERMINE);CHKERRQ(ierr); 62 ierr = VecSetBlockSize(*g,dd->w);CHKERRQ(ierr); 63 ierr = VecSetType(*g,da->vectype);CHKERRQ(ierr); 64 ierr = VecSetDM(*g, da);CHKERRQ(ierr); 65 #if defined(PETSC_HAVE_MATLAB_ENGINE) 66 if (dd->w == 1 && dd->dim == 2) { 67 ierr = PetscObjectComposeFunction((PetscObject)*g,"PetscMatlabEnginePut_C",VecMatlabEnginePut_DA2d);CHKERRQ(ierr); 68 } 69 #endif 70 } 71 PetscFunctionReturn(0); 72 } 73 74 #undef __FUNCT__ 75 #define __FUNCT__ "DMDAGetNumCells" 76 PetscErrorCode DMDAGetNumCells(DM dm, PetscInt *numCellsX, PetscInt *numCellsY, PetscInt *numCellsZ, PetscInt *numCells) 77 { 78 DM_DA *da = (DM_DA*) dm->data; 79 const PetscInt dim = da->dim; 80 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 81 const PetscInt nC = (mx)*(dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1); 82 83 PetscFunctionBegin; 84 if (numCellsX) { 85 PetscValidIntPointer(numCellsX,2); 86 *numCellsX = mx; 87 } 88 if (numCellsY) { 89 PetscValidIntPointer(numCellsX,3); 90 *numCellsY = my; 91 } 92 if (numCellsZ) { 93 PetscValidIntPointer(numCellsX,4); 94 *numCellsZ = mz; 95 } 96 if (numCells) { 97 PetscValidIntPointer(numCells,5); 98 *numCells = nC; 99 } 100 PetscFunctionReturn(0); 101 } 102 103 #undef __FUNCT__ 104 #define __FUNCT__ "DMDAGetNumVertices" 105 PetscErrorCode DMDAGetNumVertices(DM dm, PetscInt *numVerticesX, PetscInt *numVerticesY, PetscInt *numVerticesZ, PetscInt *numVertices) 106 { 107 DM_DA *da = (DM_DA*) dm->data; 108 const PetscInt dim = da->dim; 109 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 110 const PetscInt nVx = mx+1; 111 const PetscInt nVy = dim > 1 ? (my+1) : 1; 112 const PetscInt nVz = dim > 2 ? (mz+1) : 1; 113 const PetscInt nV = nVx*nVy*nVz; 114 115 PetscFunctionBegin; 116 if (numVerticesX) { 117 PetscValidIntPointer(numVerticesX,2); 118 *numVerticesX = nVx; 119 } 120 if (numVerticesY) { 121 PetscValidIntPointer(numVerticesY,3); 122 *numVerticesY = nVy; 123 } 124 if (numVerticesZ) { 125 PetscValidIntPointer(numVerticesZ,4); 126 *numVerticesZ = nVz; 127 } 128 if (numVertices) { 129 PetscValidIntPointer(numVertices,5); 130 *numVertices = nV; 131 } 132 PetscFunctionReturn(0); 133 } 134 135 #undef __FUNCT__ 136 #define __FUNCT__ "DMDAGetNumFaces" 137 PetscErrorCode DMDAGetNumFaces(DM dm, PetscInt *numXFacesX, PetscInt *numXFaces, PetscInt *numYFacesY, PetscInt *numYFaces, PetscInt *numZFacesZ, PetscInt *numZFaces) 138 { 139 DM_DA *da = (DM_DA*) dm->data; 140 const PetscInt dim = da->dim; 141 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 142 const PetscInt nxF = (dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1); 143 const PetscInt nXF = (mx+1)*nxF; 144 const PetscInt nyF = mx*(dim > 2 ? mz : 1); 145 const PetscInt nYF = dim > 1 ? (my+1)*nyF : 0; 146 const PetscInt nzF = mx*(dim > 1 ? my : 0); 147 const PetscInt nZF = dim > 2 ? (mz+1)*nzF : 0; 148 149 PetscFunctionBegin; 150 if (numXFacesX) { 151 PetscValidIntPointer(numXFacesX,2); 152 *numXFacesX = nxF; 153 } 154 if (numXFaces) { 155 PetscValidIntPointer(numXFaces,3); 156 *numXFaces = nXF; 157 } 158 if (numYFacesY) { 159 PetscValidIntPointer(numYFacesY,4); 160 *numYFacesY = nyF; 161 } 162 if (numYFaces) { 163 PetscValidIntPointer(numYFaces,5); 164 *numYFaces = nYF; 165 } 166 if (numZFacesZ) { 167 PetscValidIntPointer(numZFacesZ,6); 168 *numZFacesZ = nzF; 169 } 170 if (numZFaces) { 171 PetscValidIntPointer(numZFaces,7); 172 *numZFaces = nZF; 173 } 174 PetscFunctionReturn(0); 175 } 176 177 #undef __FUNCT__ 178 #define __FUNCT__ "DMDAGetHeightStratum" 179 PetscErrorCode DMDAGetHeightStratum(DM dm, PetscInt height, PetscInt *pStart, PetscInt *pEnd) 180 { 181 DM_DA *da = (DM_DA*) dm->data; 182 const PetscInt dim = da->dim; 183 PetscInt nC, nV, nXF, nYF, nZF; 184 PetscErrorCode ierr; 185 186 PetscFunctionBegin; 187 if (pStart) PetscValidIntPointer(pStart,3); 188 if (pEnd) PetscValidIntPointer(pEnd,4); 189 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); 190 ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr); 191 ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr); 192 if (height == 0) { 193 /* Cells */ 194 if (pStart) *pStart = 0; 195 if (pEnd) *pEnd = nC; 196 } else if (height == 1) { 197 /* Faces */ 198 if (pStart) *pStart = nC+nV; 199 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 200 } else if (height == dim) { 201 /* Vertices */ 202 if (pStart) *pStart = nC; 203 if (pEnd) *pEnd = nC+nV; 204 } else if (height < 0) { 205 /* All points */ 206 if (pStart) *pStart = 0; 207 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 208 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of height %d in the DA", height); 209 PetscFunctionReturn(0); 210 } 211 212 #undef __FUNCT__ 213 #define __FUNCT__ "DMDAGetDepthStratum" 214 PetscErrorCode DMDAGetDepthStratum(DM dm, PetscInt depth, PetscInt *pStart, PetscInt *pEnd) 215 { 216 DM_DA *da = (DM_DA*) dm->data; 217 const PetscInt dim = da->dim; 218 PetscInt nC, nV, nXF, nYF, nZF; 219 PetscErrorCode ierr; 220 221 PetscFunctionBegin; 222 if (pStart) PetscValidIntPointer(pStart,3); 223 if (pEnd) PetscValidIntPointer(pEnd,4); 224 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); 225 ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr); 226 ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr); 227 if (depth == dim) { 228 /* Cells */ 229 if (pStart) *pStart = 0; 230 if (pEnd) *pEnd = nC; 231 } else if (depth == dim-1) { 232 /* Faces */ 233 if (pStart) *pStart = nC+nV; 234 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 235 } else if (depth == 0) { 236 /* Vertices */ 237 if (pStart) *pStart = nC; 238 if (pEnd) *pEnd = nC+nV; 239 } else if (depth < 0) { 240 /* All points */ 241 if (pStart) *pStart = 0; 242 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 243 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of depth %d in the DA", depth); 244 PetscFunctionReturn(0); 245 } 246 247 #undef __FUNCT__ 248 #define __FUNCT__ "DMDAGetConeSize" 249 PetscErrorCode DMDAGetConeSize(DM dm, PetscInt p, PetscInt *coneSize) 250 { 251 DM_DA *da = (DM_DA*) dm->data; 252 const PetscInt dim = da->dim; 253 PetscInt nC, nV, nXF, nYF, nZF; 254 PetscErrorCode ierr; 255 256 PetscFunctionBegin; 257 *coneSize = 0; 258 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); 259 ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr); 260 ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr); 261 switch (dim) { 262 case 2: 263 if (p >= 0) { 264 if (p < nC) { 265 *coneSize = 4; 266 } else if (p < nC+nV) { 267 *coneSize = 0; 268 } else if (p < nC+nV+nXF+nYF+nZF) { 269 *coneSize = 2; 270 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF); 271 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p); 272 break; 273 case 3: 274 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D"); 275 break; 276 } 277 PetscFunctionReturn(0); 278 } 279 280 #undef __FUNCT__ 281 #define __FUNCT__ "DMDAGetCone" 282 PetscErrorCode DMDAGetCone(DM dm, PetscInt p, PetscInt *cone[]) 283 { 284 DM_DA *da = (DM_DA*) dm->data; 285 const PetscInt dim = da->dim; 286 PetscInt nCx, nCy, nCz, nC, nVx, nVy, nVz, nV, nxF, nyF, nzF, nXF, nYF, nZF; 287 PetscErrorCode ierr; 288 289 PetscFunctionBegin; 290 if (!cone) {ierr = DMGetWorkArray(dm, 6, PETSC_INT, cone);CHKERRQ(ierr);} 291 ierr = DMDAGetNumCells(dm, &nCx, &nCy, &nCz, &nC);CHKERRQ(ierr); 292 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr); 293 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 294 switch (dim) { 295 case 2: 296 if (p >= 0) { 297 if (p < nC) { 298 const PetscInt cy = p / nCx; 299 const PetscInt cx = p % nCx; 300 301 (*cone)[0] = cy*nxF + cx + nC+nV; 302 (*cone)[1] = cx*nyF + cy + nyF + nC+nV+nXF; 303 (*cone)[2] = cy*nxF + cx + nxF + nC+nV; 304 (*cone)[3] = cx*nyF + cy + nC+nV+nXF; 305 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do cell cones"); 306 } else if (p < nC+nV) { 307 } else if (p < nC+nV+nXF) { 308 const PetscInt fy = (p - nC+nV) / nxF; 309 const PetscInt fx = (p - nC+nV) % nxF; 310 311 (*cone)[0] = fy*nVx + fx + nC; 312 (*cone)[1] = fy*nVx + fx + 1 + nC; 313 } else if (p < nC+nV+nXF+nYF) { 314 const PetscInt fx = (p - nC+nV+nXF) / nyF; 315 const PetscInt fy = (p - nC+nV+nXF) % nyF; 316 317 (*cone)[0] = fy*nVx + fx + nC; 318 (*cone)[1] = fy*nVx + fx + nVx + nC; 319 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF); 320 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p); 321 break; 322 case 3: 323 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D"); 324 break; 325 } 326 PetscFunctionReturn(0); 327 } 328 329 #undef __FUNCT__ 330 #define __FUNCT__ "DMDARestoreCone" 331 PetscErrorCode DMDARestoreCone(DM dm, PetscInt p, PetscInt *cone[]) 332 { 333 PetscErrorCode ierr; 334 335 PetscFunctionBegin; 336 ierr = DMGetWorkArray(dm, 6, PETSC_INT, cone);CHKERRQ(ierr); 337 PetscFunctionReturn(0); 338 } 339 340 #undef __FUNCT__ 341 #define __FUNCT__ "DMDACreateSection" 342 /*@C 343 DMDACreateSection - Create a PetscSection inside the DMDA that describes data layout. This allows multiple fields with 344 different numbers of dofs on vertices, cells, and faces in each direction. 345 346 Input Parameters: 347 + dm- The DMDA 348 . numFields - The number of fields 349 . numComp - The number of components in each field 350 . numDof - The number of dofs per dimension for each field 351 . numFaceDof - The number of dofs per face for each field and direction, or NULL 352 353 Level: developer 354 355 Note: 356 The default DMDA numbering is as follows: 357 358 - Cells: [0, nC) 359 - Vertices: [nC, nC+nV) 360 - X-Faces: [nC+nV, nC+nV+nXF) normal is +- x-dir 361 - Y-Faces: [nC+nV+nXF, nC+nV+nXF+nYF) normal is +- y-dir 362 - Z-Faces: [nC+nV+nXF+nYF, nC+nV+nXF+nYF+nZF) normal is +- z-dir 363 364 We interpret the default DMDA partition as a cell partition, and the data assignment as a cell assignment. 365 @*/ 366 PetscErrorCode DMDACreateSection(DM dm, PetscInt numComp[], PetscInt numDof[], PetscInt numFaceDof[], PetscSection *s) 367 { 368 DM_DA *da = (DM_DA*) dm->data; 369 PetscSection section; 370 const PetscInt dim = da->dim; 371 PetscInt numFields, numVertexTotDof = 0, numCellTotDof = 0, numFaceTotDof[3] = {0, 0, 0}; 372 PetscSF sf; 373 PetscMPIInt rank; 374 const PetscMPIInt *neighbors; 375 PetscInt *localPoints; 376 PetscSFNode *remotePoints; 377 PetscInt nleaves = 0, nleavesCheck = 0, nL = 0; 378 PetscInt nC, nVx, nVy, nVz, nV, nxF, nXF, nyF, nYF, nzF, nZF; 379 PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart, zfEnd; 380 PetscInt f, v, c, xf, yf, zf, xn, yn, zn; 381 PetscErrorCode ierr; 382 383 PetscFunctionBegin; 384 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 385 PetscValidPointer(s, 4); 386 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 387 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); 388 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr); 389 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 390 ierr = DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);CHKERRQ(ierr); 391 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 392 ierr = DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 393 ierr = DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);CHKERRQ(ierr); 394 xfStart = vEnd; xfEnd = xfStart+nXF; 395 yfStart = xfEnd; yfEnd = yfStart+nYF; 396 zfStart = yfEnd; zfEnd = zfStart+nZF; 397 if (zfEnd != fEnd) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid face end %d, should be %d", zfEnd, fEnd); 398 /* Create local section */ 399 ierr = DMDAGetInfo(dm, 0,0,0,0,0,0,0, &numFields, 0,0,0,0,0);CHKERRQ(ierr); 400 for (f = 0; f < numFields; ++f) { 401 numVertexTotDof += numDof[f*(dim+1)+0]; 402 numCellTotDof += numDof[f*(dim+1)+dim]; 403 numFaceTotDof[0] += dim > 0 ? (numFaceDof ? numFaceDof[f*dim+0] : numDof[f*(dim+1)+dim-1]) : 0; 404 numFaceTotDof[1] += dim > 1 ? (numFaceDof ? numFaceDof[f*dim+1] : numDof[f*(dim+1)+dim-1]) : 0; 405 numFaceTotDof[2] += dim > 2 ? (numFaceDof ? numFaceDof[f*dim+2] : numDof[f*(dim+1)+dim-1]) : 0; 406 } 407 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion);CHKERRQ(ierr); 408 if (numFields > 0) { 409 ierr = PetscSectionSetNumFields(section, numFields);CHKERRQ(ierr); 410 for (f = 0; f < numFields; ++f) { 411 const char *name; 412 413 ierr = DMDAGetFieldName(dm, f, &name);CHKERRQ(ierr); 414 ierr = PetscSectionSetFieldName(section, f, name ? name : "Field");CHKERRQ(ierr); 415 if (numComp) { 416 ierr = PetscSectionSetFieldComponents(section, f, numComp[f]);CHKERRQ(ierr); 417 } 418 } 419 } 420 ierr = PetscSectionSetChart(section, pStart, pEnd);CHKERRQ(ierr); 421 for (v = vStart; v < vEnd; ++v) { 422 for (f = 0; f < numFields; ++f) { 423 ierr = PetscSectionSetFieldDof(section, v, f, numDof[f*(dim+1)+0]);CHKERRQ(ierr); 424 } 425 ierr = PetscSectionSetDof(section, v, numVertexTotDof);CHKERRQ(ierr); 426 } 427 for (xf = xfStart; xf < xfEnd; ++xf) { 428 for (f = 0; f < numFields; ++f) { 429 ierr = PetscSectionSetFieldDof(section, xf, f, numFaceDof ? numFaceDof[f*dim+0] : numDof[f*(dim+1)+dim-1]);CHKERRQ(ierr); 430 } 431 ierr = PetscSectionSetDof(section, xf, numFaceTotDof[0]);CHKERRQ(ierr); 432 } 433 for (yf = yfStart; yf < yfEnd; ++yf) { 434 for (f = 0; f < numFields; ++f) { 435 ierr = PetscSectionSetFieldDof(section, yf, f, numFaceDof ? numFaceDof[f*dim+1] : numDof[f*(dim+1)+dim-1]);CHKERRQ(ierr); 436 } 437 ierr = PetscSectionSetDof(section, yf, numFaceTotDof[1]);CHKERRQ(ierr); 438 } 439 for (zf = zfStart; zf < zfEnd; ++zf) { 440 for (f = 0; f < numFields; ++f) { 441 ierr = PetscSectionSetFieldDof(section, zf, f, numFaceDof ? numFaceDof[f*dim+2] : numDof[f*(dim+1)+dim-1]);CHKERRQ(ierr); 442 } 443 ierr = PetscSectionSetDof(section, zf, numFaceTotDof[2]);CHKERRQ(ierr); 444 } 445 for (c = cStart; c < cEnd; ++c) { 446 for (f = 0; f < numFields; ++f) { 447 ierr = PetscSectionSetFieldDof(section, c, f, numDof[f*(dim+1)+dim]);CHKERRQ(ierr); 448 } 449 ierr = PetscSectionSetDof(section, c, numCellTotDof);CHKERRQ(ierr); 450 } 451 ierr = PetscSectionSetUp(section);CHKERRQ(ierr); 452 /* Create mesh point SF */ 453 ierr = DMDAGetNeighbors(dm, &neighbors);CHKERRQ(ierr); 454 for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) { 455 for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) { 456 for (xn = 0; xn < 3; ++xn) { 457 const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0; 458 const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn]; 459 460 if (neighbor >= 0 && neighbor < rank) { 461 nleaves += (!xp ? nVx : 1) * (!yp ? nVy : 1) * (!zp ? nVz : 1); /* vertices */ 462 if (xp && !yp && !zp) { 463 nleaves += nxF; /* x faces */ 464 } else if (yp && !zp && !xp) { 465 nleaves += nyF; /* y faces */ 466 } else if (zp && !xp && !yp) { 467 nleaves += nzF; /* z faces */ 468 } 469 } 470 } 471 } 472 } 473 ierr = PetscMalloc2(nleaves,&localPoints,nleaves,&remotePoints);CHKERRQ(ierr); 474 for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) { 475 for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) { 476 for (xn = 0; xn < 3; ++xn) { 477 const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0; 478 const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn]; 479 PetscInt xv, yv, zv; 480 481 if (neighbor >= 0 && neighbor < rank) { 482 if (xp < 0) { /* left */ 483 if (yp < 0) { /* bottom */ 484 if (zp < 0) { /* back */ 485 const PetscInt localVertex = ( 0*nVy + 0)*nVx + 0 + nC; 486 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 487 nleavesCheck += 1; /* left bottom back vertex */ 488 489 localPoints[nL] = localVertex; 490 remotePoints[nL].rank = neighbor; 491 remotePoints[nL].index = remoteVertex; 492 ++nL; 493 } else if (zp > 0) { /* front */ 494 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; 495 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 496 nleavesCheck += 1; /* left bottom front vertex */ 497 498 localPoints[nL] = localVertex; 499 remotePoints[nL].rank = neighbor; 500 remotePoints[nL].index = remoteVertex; 501 ++nL; 502 } else { 503 nleavesCheck += nVz; /* left bottom vertices */ 504 for (zv = 0; zv < nVz; ++zv, ++nL) { 505 const PetscInt localVertex = (zv*nVy + 0)*nVx + 0 + nC; 506 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 507 508 localPoints[nL] = localVertex; 509 remotePoints[nL].rank = neighbor; 510 remotePoints[nL].index = remoteVertex; 511 } 512 } 513 } else if (yp > 0) { /* top */ 514 if (zp < 0) { /* back */ 515 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; 516 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 517 nleavesCheck += 1; /* left top back vertex */ 518 519 localPoints[nL] = localVertex; 520 remotePoints[nL].rank = neighbor; 521 remotePoints[nL].index = remoteVertex; 522 ++nL; 523 } else if (zp > 0) { /* front */ 524 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; 525 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 526 nleavesCheck += 1; /* left top front vertex */ 527 528 localPoints[nL] = localVertex; 529 remotePoints[nL].rank = neighbor; 530 remotePoints[nL].index = remoteVertex; 531 ++nL; 532 } else { 533 nleavesCheck += nVz; /* left top vertices */ 534 for (zv = 0; zv < nVz; ++zv, ++nL) { 535 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; 536 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 537 538 localPoints[nL] = localVertex; 539 remotePoints[nL].rank = neighbor; 540 remotePoints[nL].index = remoteVertex; 541 } 542 } 543 } else { 544 if (zp < 0) { /* back */ 545 nleavesCheck += nVy; /* left back vertices */ 546 for (yv = 0; yv < nVy; ++yv, ++nL) { 547 const PetscInt localVertex = ( 0*nVy + yv)*nVx + 0 + nC; 548 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 549 550 localPoints[nL] = localVertex; 551 remotePoints[nL].rank = neighbor; 552 remotePoints[nL].index = remoteVertex; 553 } 554 } else if (zp > 0) { /* front */ 555 nleavesCheck += nVy; /* left front vertices */ 556 for (yv = 0; yv < nVy; ++yv, ++nL) { 557 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; 558 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 559 560 localPoints[nL] = localVertex; 561 remotePoints[nL].rank = neighbor; 562 remotePoints[nL].index = remoteVertex; 563 } 564 } else { 565 nleavesCheck += nVy*nVz; /* left vertices */ 566 for (zv = 0; zv < nVz; ++zv) { 567 for (yv = 0; yv < nVy; ++yv, ++nL) { 568 const PetscInt localVertex = (zv*nVy + yv)*nVx + 0 + nC; 569 const PetscInt remoteVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 570 571 localPoints[nL] = localVertex; 572 remotePoints[nL].rank = neighbor; 573 remotePoints[nL].index = remoteVertex; 574 } 575 } 576 nleavesCheck += nxF; /* left faces */ 577 for (xf = 0; xf < nxF; ++xf, ++nL) { 578 /* THIS IS WRONG */ 579 const PetscInt localFace = 0 + nC+nV; 580 const PetscInt remoteFace = 0 + nC+nV; 581 582 localPoints[nL] = localFace; 583 remotePoints[nL].rank = neighbor; 584 remotePoints[nL].index = remoteFace; 585 } 586 } 587 } 588 } else if (xp > 0) { /* right */ 589 if (yp < 0) { /* bottom */ 590 if (zp < 0) { /* back */ 591 const PetscInt localVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; 592 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 593 nleavesCheck += 1; /* right bottom back vertex */ 594 595 localPoints[nL] = localVertex; 596 remotePoints[nL].rank = neighbor; 597 remotePoints[nL].index = remoteVertex; 598 ++nL; 599 } else if (zp > 0) { /* front */ 600 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; 601 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 602 nleavesCheck += 1; /* right bottom front vertex */ 603 604 localPoints[nL] = localVertex; 605 remotePoints[nL].rank = neighbor; 606 remotePoints[nL].index = remoteVertex; 607 ++nL; 608 } else { 609 nleavesCheck += nVz; /* right bottom vertices */ 610 for (zv = 0; zv < nVz; ++zv, ++nL) { 611 const PetscInt localVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; 612 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 613 614 localPoints[nL] = localVertex; 615 remotePoints[nL].rank = neighbor; 616 remotePoints[nL].index = remoteVertex; 617 } 618 } 619 } else if (yp > 0) { /* top */ 620 if (zp < 0) { /* back */ 621 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; 622 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 623 nleavesCheck += 1; /* right top back vertex */ 624 625 localPoints[nL] = localVertex; 626 remotePoints[nL].rank = neighbor; 627 remotePoints[nL].index = remoteVertex; 628 ++nL; 629 } else if (zp > 0) { /* front */ 630 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; 631 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 632 nleavesCheck += 1; /* right top front vertex */ 633 634 localPoints[nL] = localVertex; 635 remotePoints[nL].rank = neighbor; 636 remotePoints[nL].index = remoteVertex; 637 ++nL; 638 } else { 639 nleavesCheck += nVz; /* right top vertices */ 640 for (zv = 0; zv < nVz; ++zv, ++nL) { 641 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; 642 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 643 644 localPoints[nL] = localVertex; 645 remotePoints[nL].rank = neighbor; 646 remotePoints[nL].index = remoteVertex; 647 } 648 } 649 } else { 650 if (zp < 0) { /* back */ 651 nleavesCheck += nVy; /* right back vertices */ 652 for (yv = 0; yv < nVy; ++yv, ++nL) { 653 const PetscInt localVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; 654 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 655 656 localPoints[nL] = localVertex; 657 remotePoints[nL].rank = neighbor; 658 remotePoints[nL].index = remoteVertex; 659 } 660 } else if (zp > 0) { /* front */ 661 nleavesCheck += nVy; /* right front vertices */ 662 for (yv = 0; yv < nVy; ++yv, ++nL) { 663 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; 664 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 665 666 localPoints[nL] = localVertex; 667 remotePoints[nL].rank = neighbor; 668 remotePoints[nL].index = remoteVertex; 669 } 670 } else { 671 nleavesCheck += nVy*nVz; /* right vertices */ 672 for (zv = 0; zv < nVz; ++zv) { 673 for (yv = 0; yv < nVy; ++yv, ++nL) { 674 const PetscInt localVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; 675 const PetscInt remoteVertex = (zv*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 676 677 localPoints[nL] = localVertex; 678 remotePoints[nL].rank = neighbor; 679 remotePoints[nL].index = remoteVertex; 680 } 681 } 682 nleavesCheck += nxF; /* right faces */ 683 for (xf = 0; xf < nxF; ++xf, ++nL) { 684 /* THIS IS WRONG */ 685 const PetscInt localFace = 0 + nC+nV; 686 const PetscInt remoteFace = 0 + nC+nV; 687 688 localPoints[nL] = localFace; 689 remotePoints[nL].rank = neighbor; 690 remotePoints[nL].index = remoteFace; 691 } 692 } 693 } 694 } else { 695 if (yp < 0) { /* bottom */ 696 if (zp < 0) { /* back */ 697 nleavesCheck += nVx; /* bottom back vertices */ 698 for (xv = 0; xv < nVx; ++xv, ++nL) { 699 const PetscInt localVertex = ( 0*nVy + 0)*nVx + xv + nC; 700 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 701 702 localPoints[nL] = localVertex; 703 remotePoints[nL].rank = neighbor; 704 remotePoints[nL].index = remoteVertex; 705 } 706 } else if (zp > 0) { /* front */ 707 nleavesCheck += nVx; /* bottom front vertices */ 708 for (xv = 0; xv < nVx; ++xv, ++nL) { 709 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; 710 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 711 712 localPoints[nL] = localVertex; 713 remotePoints[nL].rank = neighbor; 714 remotePoints[nL].index = remoteVertex; 715 } 716 } else { 717 nleavesCheck += nVx*nVz; /* bottom vertices */ 718 for (zv = 0; zv < nVz; ++zv) { 719 for (xv = 0; xv < nVx; ++xv, ++nL) { 720 const PetscInt localVertex = (zv*nVy + 0)*nVx + xv + nC; 721 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 722 723 localPoints[nL] = localVertex; 724 remotePoints[nL].rank = neighbor; 725 remotePoints[nL].index = remoteVertex; 726 } 727 } 728 nleavesCheck += nyF; /* bottom faces */ 729 for (yf = 0; yf < nyF; ++yf, ++nL) { 730 /* THIS IS WRONG */ 731 const PetscInt localFace = 0 + nC+nV; 732 const PetscInt remoteFace = 0 + nC+nV; 733 734 localPoints[nL] = localFace; 735 remotePoints[nL].rank = neighbor; 736 remotePoints[nL].index = remoteFace; 737 } 738 } 739 } else if (yp > 0) { /* top */ 740 if (zp < 0) { /* back */ 741 nleavesCheck += nVx; /* top back vertices */ 742 for (xv = 0; xv < nVx; ++xv, ++nL) { 743 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; 744 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 745 746 localPoints[nL] = localVertex; 747 remotePoints[nL].rank = neighbor; 748 remotePoints[nL].index = remoteVertex; 749 } 750 } else if (zp > 0) { /* front */ 751 nleavesCheck += nVx; /* top front vertices */ 752 for (xv = 0; xv < nVx; ++xv, ++nL) { 753 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; 754 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 755 756 localPoints[nL] = localVertex; 757 remotePoints[nL].rank = neighbor; 758 remotePoints[nL].index = remoteVertex; 759 } 760 } else { 761 nleavesCheck += nVx*nVz; /* top vertices */ 762 for (zv = 0; zv < nVz; ++zv) { 763 for (xv = 0; xv < nVx; ++xv, ++nL) { 764 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + xv + nC; 765 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 766 767 localPoints[nL] = localVertex; 768 remotePoints[nL].rank = neighbor; 769 remotePoints[nL].index = remoteVertex; 770 } 771 } 772 nleavesCheck += nyF; /* top faces */ 773 for (yf = 0; yf < nyF; ++yf, ++nL) { 774 /* THIS IS WRONG */ 775 const PetscInt localFace = 0 + nC+nV; 776 const PetscInt remoteFace = 0 + nC+nV; 777 778 localPoints[nL] = localFace; 779 remotePoints[nL].rank = neighbor; 780 remotePoints[nL].index = remoteFace; 781 } 782 } 783 } else { 784 if (zp < 0) { /* back */ 785 nleavesCheck += nVx*nVy; /* back vertices */ 786 for (yv = 0; yv < nVy; ++yv) { 787 for (xv = 0; xv < nVx; ++xv, ++nL) { 788 const PetscInt localVertex = ( 0*nVy + yv)*nVx + xv + nC; 789 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 790 791 localPoints[nL] = localVertex; 792 remotePoints[nL].rank = neighbor; 793 remotePoints[nL].index = remoteVertex; 794 } 795 } 796 nleavesCheck += nzF; /* back faces */ 797 for (zf = 0; zf < nzF; ++zf, ++nL) { 798 /* THIS IS WRONG */ 799 const PetscInt localFace = 0 + nC+nV; 800 const PetscInt remoteFace = 0 + nC+nV; 801 802 localPoints[nL] = localFace; 803 remotePoints[nL].rank = neighbor; 804 remotePoints[nL].index = remoteFace; 805 } 806 } else if (zp > 0) { /* front */ 807 nleavesCheck += nVx*nVy; /* front vertices */ 808 for (yv = 0; yv < nVy; ++yv) { 809 for (xv = 0; xv < nVx; ++xv, ++nL) { 810 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; 811 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 812 813 localPoints[nL] = localVertex; 814 remotePoints[nL].rank = neighbor; 815 remotePoints[nL].index = remoteVertex; 816 } 817 } 818 nleavesCheck += nzF; /* front faces */ 819 for (zf = 0; zf < nzF; ++zf, ++nL) { 820 /* THIS IS WRONG */ 821 const PetscInt localFace = 0 + nC+nV; 822 const PetscInt remoteFace = 0 + nC+nV; 823 824 localPoints[nL] = localFace; 825 remotePoints[nL].rank = neighbor; 826 remotePoints[nL].index = remoteFace; 827 } 828 } else { 829 /* Nothing is shared from the interior */ 830 } 831 } 832 } 833 } 834 } 835 } 836 } 837 /* TODO: Remove duplication in leaf determination */ 838 if (nleaves != nleavesCheck) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "The number of leaves %d did not match the number of remote leaves %d", nleaves, nleavesCheck); 839 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf);CHKERRQ(ierr); 840 ierr = PetscSFSetGraph(sf, pEnd, nleaves, localPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 841 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 842 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 843 *s = section; 844 PetscFunctionReturn(0); 845 } 846 847 #undef __FUNCT__ 848 #define __FUNCT__ "DMDASetVertexCoordinates" 849 PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu) 850 { 851 DM_DA *da = (DM_DA *) dm->data; 852 Vec coordinates; 853 PetscSection section; 854 PetscScalar *coords; 855 PetscReal h[3]; 856 PetscInt dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k; 857 PetscErrorCode ierr; 858 859 PetscFunctionBegin; 860 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 861 ierr = DMDAGetInfo(dm, &dim, &M, &N, &P, 0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 862 h[0] = (xu - xl)/M; 863 h[1] = (yu - yl)/N; 864 h[2] = (zu - zl)/P; 865 ierr = DMDAGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 866 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr); 867 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ion);CHKERRQ(ierr); 868 ierr = PetscSectionSetNumFields(section, 1);CHKERRQ(ierr); 869 ierr = PetscSectionSetFieldComponents(section, 0, dim);CHKERRQ(ierr); 870 ierr = PetscSectionSetChart(section, vStart, vEnd);CHKERRQ(ierr); 871 for (v = vStart; v < vEnd; ++v) { 872 ierr = PetscSectionSetDof(section, v, dim);CHKERRQ(ierr); 873 } 874 ierr = PetscSectionSetUp(section);CHKERRQ(ierr); 875 ierr = PetscSectionGetStorageSize(section, &size);CHKERRQ(ierr); 876 ierr = VecCreateSeq(PETSC_COMM_SELF, size, &coordinates);CHKERRQ(ierr); 877 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 878 for (k = 0; k < nVz; ++k) { 879 PetscInt ind[3] = {0, 0, k + da->zs}, d, off; 880 881 for (j = 0; j < nVy; ++j) { 882 ind[1] = j + da->ys; 883 for (i = 0; i < nVx; ++i) { 884 const PetscInt vertex = (k*nVy + j)*nVx + i + vStart; 885 886 ierr = PetscSectionGetOffset(section, vertex, &off);CHKERRQ(ierr); 887 ind[0] = i + da->xs; 888 for (d = 0; d < dim; ++d) { 889 coords[off+d] = h[d]*ind[d]; 890 } 891 } 892 } 893 } 894 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 895 ierr = DMSetCoordinateSection(dm, section);CHKERRQ(ierr); 896 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 897 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr); 898 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 899 PetscFunctionReturn(0); 900 } 901 PetscFunctionReturn(0); 902 } 903 904 /* ------------------------------------------------------------------- */ 905 906 #undef __FUNCT__ 907 #define __FUNCT__ "DMDAGetArray" 908 /*@C 909 DMDAGetArray - Gets a work array for a DMDA 910 911 Input Parameter: 912 + da - information about my local patch 913 - ghosted - do you want arrays for the ghosted or nonghosted patch 914 915 Output Parameters: 916 . vptr - array data structured 917 918 Note: The vector values are NOT initialized and may have garbage in them, so you may need 919 to zero them. 920 921 Level: advanced 922 923 .seealso: DMDARestoreArray() 924 925 @*/ 926 PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr) 927 { 928 PetscErrorCode ierr; 929 PetscInt j,i,xs,ys,xm,ym,zs,zm; 930 char *iarray_start; 931 void **iptr = (void**)vptr; 932 DM_DA *dd = (DM_DA*)da->data; 933 934 PetscFunctionBegin; 935 PetscValidHeaderSpecific(da,DM_CLASSID,1); 936 if (ghosted) { 937 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 938 if (dd->arrayghostedin[i]) { 939 *iptr = dd->arrayghostedin[i]; 940 iarray_start = (char*)dd->startghostedin[i]; 941 dd->arrayghostedin[i] = NULL; 942 dd->startghostedin[i] = NULL; 943 944 goto done; 945 } 946 } 947 xs = dd->Xs; 948 ys = dd->Ys; 949 zs = dd->Zs; 950 xm = dd->Xe-dd->Xs; 951 ym = dd->Ye-dd->Ys; 952 zm = dd->Ze-dd->Zs; 953 } else { 954 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 955 if (dd->arrayin[i]) { 956 *iptr = dd->arrayin[i]; 957 iarray_start = (char*)dd->startin[i]; 958 dd->arrayin[i] = NULL; 959 dd->startin[i] = NULL; 960 961 goto done; 962 } 963 } 964 xs = dd->xs; 965 ys = dd->ys; 966 zs = dd->zs; 967 xm = dd->xe-dd->xs; 968 ym = dd->ye-dd->ys; 969 zm = dd->ze-dd->zs; 970 } 971 972 switch (dd->dim) { 973 case 1: { 974 void *ptr; 975 976 ierr = PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 977 978 ptr = (void*)(iarray_start - xs*sizeof(PetscScalar)); 979 *iptr = (void*)ptr; 980 break; 981 } 982 case 2: { 983 void **ptr; 984 985 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 986 987 ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*)); 988 for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs); 989 *iptr = (void*)ptr; 990 break; 991 } 992 case 3: { 993 void ***ptr,**bptr; 994 995 ierr = PetscMalloc((zm+1)*sizeof(void**)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 996 997 ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*)); 998 bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**)); 999 for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys); 1000 for (i=zs; i<zs+zm; i++) { 1001 for (j=ys; j<ys+ym; j++) { 1002 ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs); 1003 } 1004 } 1005 *iptr = (void*)ptr; 1006 break; 1007 } 1008 default: 1009 SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 1010 } 1011 1012 done: 1013 /* add arrays to the checked out list */ 1014 if (ghosted) { 1015 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1016 if (!dd->arrayghostedout[i]) { 1017 dd->arrayghostedout[i] = *iptr; 1018 dd->startghostedout[i] = iarray_start; 1019 break; 1020 } 1021 } 1022 } else { 1023 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1024 if (!dd->arrayout[i]) { 1025 dd->arrayout[i] = *iptr; 1026 dd->startout[i] = iarray_start; 1027 break; 1028 } 1029 } 1030 } 1031 PetscFunctionReturn(0); 1032 } 1033 1034 #undef __FUNCT__ 1035 #define __FUNCT__ "DMDARestoreArray" 1036 /*@C 1037 DMDARestoreArray - Restores an array of derivative types for a DMDA 1038 1039 Input Parameter: 1040 + da - information about my local patch 1041 . ghosted - do you want arrays for the ghosted or nonghosted patch 1042 - vptr - array data structured to be passed to ad_FormFunctionLocal() 1043 1044 Level: advanced 1045 1046 .seealso: DMDAGetArray() 1047 1048 @*/ 1049 PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr) 1050 { 1051 PetscInt i; 1052 void **iptr = (void**)vptr,*iarray_start = 0; 1053 DM_DA *dd = (DM_DA*)da->data; 1054 1055 PetscFunctionBegin; 1056 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1057 if (ghosted) { 1058 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1059 if (dd->arrayghostedout[i] == *iptr) { 1060 iarray_start = dd->startghostedout[i]; 1061 dd->arrayghostedout[i] = NULL; 1062 dd->startghostedout[i] = NULL; 1063 break; 1064 } 1065 } 1066 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1067 if (!dd->arrayghostedin[i]) { 1068 dd->arrayghostedin[i] = *iptr; 1069 dd->startghostedin[i] = iarray_start; 1070 break; 1071 } 1072 } 1073 } else { 1074 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1075 if (dd->arrayout[i] == *iptr) { 1076 iarray_start = dd->startout[i]; 1077 dd->arrayout[i] = NULL; 1078 dd->startout[i] = NULL; 1079 break; 1080 } 1081 } 1082 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1083 if (!dd->arrayin[i]) { 1084 dd->arrayin[i] = *iptr; 1085 dd->startin[i] = iarray_start; 1086 break; 1087 } 1088 } 1089 } 1090 PetscFunctionReturn(0); 1091 } 1092 1093