1 2 /* 3 Code for manipulating distributed regular arrays in parallel. 4 */ 5 6 #include <petsc-private/daimpl.h> /*I "petscdmda.h" I*/ 7 8 /* 9 This allows the DMDA vectors to properly tell MATLAB their dimensions 10 */ 11 #if defined(PETSC_HAVE_MATLAB_ENGINE) 12 #include <engine.h> /* MATLAB include file */ 13 #include <mex.h> /* MATLAB include file */ 14 EXTERN_C_BEGIN 15 #undef __FUNCT__ 16 #define __FUNCT__ "VecMatlabEnginePut_DA2d" 17 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 = PetscObjectQuery((PetscObject)vec,"DM",(PetscObject*)&da);CHKERRQ(ierr); 28 if (!da) SETERRQ(((PetscObject)vec)->comm,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 EXTERN_C_END 45 #endif 46 47 48 #undef __FUNCT__ 49 #define __FUNCT__ "DMCreateLocalVector_DA" 50 PetscErrorCode DMCreateLocalVector_DA(DM da,Vec* g) 51 { 52 PetscErrorCode ierr; 53 DM_DA *dd = (DM_DA*)da->data; 54 55 PetscFunctionBegin; 56 PetscValidHeaderSpecific(da,DM_CLASSID,1); 57 PetscValidPointer(g,2); 58 ierr = VecCreate(PETSC_COMM_SELF,g);CHKERRQ(ierr); 59 ierr = VecSetSizes(*g,dd->nlocal,PETSC_DETERMINE);CHKERRQ(ierr); 60 ierr = VecSetBlockSize(*g,dd->w);CHKERRQ(ierr); 61 ierr = VecSetType(*g,da->vectype);CHKERRQ(ierr); 62 ierr = PetscObjectCompose((PetscObject)*g,"DM",(PetscObject)da);CHKERRQ(ierr); 63 #if defined(PETSC_HAVE_MATLAB_ENGINE) 64 if (dd->w == 1 && dd->dim == 2) { 65 ierr = PetscObjectComposeFunctionDynamic((PetscObject)*g,"PetscMatlabEnginePut_C","VecMatlabEnginePut_DA2d",VecMatlabEnginePut_DA2d);CHKERRQ(ierr); 66 } 67 #endif 68 PetscFunctionReturn(0); 69 } 70 71 #undef __FUNCT__ 72 #define __FUNCT__ "DMDAGetNumCells" 73 PetscErrorCode DMDAGetNumCells(DM dm, PetscInt *numCells) 74 { 75 DM_DA *da = (DM_DA *) dm->data; 76 const PetscInt dim = da->dim; 77 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 78 const PetscInt nC = (mx )*(dim > 1 ? (my )*(dim > 2 ? (mz ) : 1) : 1); 79 80 PetscFunctionBegin; 81 if (numCells) { 82 PetscValidIntPointer(numCells,2); 83 *numCells = nC; 84 } 85 PetscFunctionReturn(0); 86 } 87 88 #undef __FUNCT__ 89 #define __FUNCT__ "DMDAGetNumVertices" 90 PetscErrorCode DMDAGetNumVertices(DM dm, PetscInt *numVerticesX, PetscInt *numVerticesY, PetscInt *numVerticesZ, PetscInt *numVertices) 91 { 92 DM_DA *da = (DM_DA *) dm->data; 93 const PetscInt dim = da->dim; 94 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 95 const PetscInt nVx = mx+1; 96 const PetscInt nVy = dim > 1 ? (my+1) : 1; 97 const PetscInt nVz = dim > 2 ? (mz+1) : 1; 98 const PetscInt nV = nVx*nVy*nVz; 99 100 PetscFunctionBegin; 101 if (numVerticesX) { 102 PetscValidIntPointer(numVerticesX,2); 103 *numVerticesX = nVx; 104 } 105 if (numVerticesY) { 106 PetscValidIntPointer(numVerticesY,3); 107 *numVerticesY = nVy; 108 } 109 if (numVerticesZ) { 110 PetscValidIntPointer(numVerticesZ,4); 111 *numVerticesZ = nVz; 112 } 113 if (numVertices) { 114 PetscValidIntPointer(numVertices,5); 115 *numVertices = nV; 116 } 117 PetscFunctionReturn(0); 118 } 119 120 #undef __FUNCT__ 121 #define __FUNCT__ "DMDAGetNumFaces" 122 PetscErrorCode DMDAGetNumFaces(DM dm, PetscInt *numXFacesX, PetscInt *numXFaces, PetscInt *numYFacesY, PetscInt *numYFaces, PetscInt *numZFacesZ, PetscInt *numZFaces) 123 { 124 DM_DA *da = (DM_DA *) dm->data; 125 const PetscInt dim = da->dim; 126 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 127 const PetscInt nxF = (dim > 1 ? (my )*(dim > 2 ? (mz ) : 1) : 1); 128 const PetscInt nXF = (mx+1)*nxF; 129 const PetscInt nyF = mx*(dim > 2 ? mz : 1); 130 const PetscInt nYF = dim > 1 ? (my+1)*nyF : 0; 131 const PetscInt nzF = mx*(dim > 1 ? my : 0); 132 const PetscInt nZF = dim > 2 ? (mz+1)*nzF : 0; 133 134 PetscFunctionBegin; 135 if (numXFacesX) { 136 PetscValidIntPointer(numXFacesX,2); 137 *numXFacesX = nxF; 138 } 139 if (numXFaces) { 140 PetscValidIntPointer(numXFaces,3); 141 *numXFaces = nXF; 142 } 143 if (numYFacesY) { 144 PetscValidIntPointer(numYFacesY,4); 145 *numYFacesY = nyF; 146 } 147 if (numYFaces) { 148 PetscValidIntPointer(numYFaces,5); 149 *numYFaces = nYF; 150 } 151 if (numZFacesZ) { 152 PetscValidIntPointer(numZFacesZ,6); 153 *numZFacesZ = nzF; 154 } 155 if (numZFaces) { 156 PetscValidIntPointer(numZFaces,7); 157 *numZFaces = nZF; 158 } 159 PetscFunctionReturn(0); 160 } 161 162 #undef __FUNCT__ 163 #define __FUNCT__ "DMDAGetHeightStratum" 164 PetscErrorCode DMDAGetHeightStratum(DM dm, PetscInt height, PetscInt *pStart, PetscInt *pEnd) 165 { 166 DM_DA *da = (DM_DA *) dm->data; 167 const PetscInt dim = da->dim; 168 PetscInt nC, nV, nXF, nYF, nZF; 169 PetscErrorCode ierr; 170 171 PetscFunctionBegin; 172 if (pStart) {PetscValidIntPointer(pStart,3);} 173 if (pEnd) {PetscValidIntPointer(pEnd,4);} 174 ierr = DMDAGetNumCells(dm, &nC);CHKERRQ(ierr); 175 ierr = DMDAGetNumVertices(dm, PETSC_NULL, PETSC_NULL, PETSC_NULL, &nV);CHKERRQ(ierr); 176 ierr = DMDAGetNumFaces(dm, PETSC_NULL, &nXF, PETSC_NULL, &nYF, PETSC_NULL, &nZF);CHKERRQ(ierr); 177 if (height == 0) { 178 /* Cells */ 179 if (pStart) {*pStart = 0;} 180 if (pEnd) {*pEnd = nC;} 181 } else if (height == 1) { 182 /* Faces */ 183 if (pStart) {*pStart = nC+nV;} 184 if (pEnd) {*pEnd = nC+nV+nXF+nYF+nZF;} 185 } else if (height == dim) { 186 /* Vertices */ 187 if (pStart) {*pStart = nC;} 188 if (pEnd) {*pEnd = nC+nV;} 189 } else if (height < 0) { 190 /* All points */ 191 if (pStart) {*pStart = 0;} 192 if (pEnd) {*pEnd = nC+nV+nXF+nYF+nZF;} 193 } else { 194 SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "No points of height %d in the DA", height); 195 } 196 PetscFunctionReturn(0); 197 } 198 199 #undef __FUNCT__ 200 #define __FUNCT__ "DMDACreateSection" 201 /*@C 202 DMDACreateSection - Create a PetscSection inside the DMDA that describes data layout. This allows multiple fields with 203 different numbers of dofs on vertices, cells, and faces in each direction. 204 205 Input Parameters: 206 + dm- The DMDA 207 . numFields - The number of fields 208 . numComp - The number of components in each field, or PETSC_NULL for 1 209 . numVertexDof - The number of dofs per vertex for each field, or PETSC_NULL 210 . numFaceDof - The number of dofs per face for each field and direction, or PETSC_NULL 211 - numCellDof - The number of dofs per cell for each field, or PETSC_NULL 212 213 Level: developer 214 215 Note: 216 The default DMDA numbering is as follows: 217 218 - Cells: [0, nC) 219 - Vertices: [nC, nC+nV) 220 - X-Faces: [nC+nV, nC+nV+nXF) normal is +- x-dir 221 - Y-Faces: [nC+nV+nXF, nC+nV+nXF+nYF) normal is +- y-dir 222 - Z-Faces: [nC+nV+nXF+nYF, nC+nV+nXF+nYF+nZF) normal is +- z-dir 223 224 We interpret the default DMDA partition as a cell partition, and the data assignment as a cell assignment. 225 @*/ 226 PetscErrorCode DMDACreateSection(DM dm, PetscInt numComp[], PetscInt numVertexDof[], PetscInt numFaceDof[], PetscInt numCellDof[]) 227 { 228 DM_DA *da = (DM_DA *) dm->data; 229 const PetscInt dim = da->dim; 230 PetscInt numFields, numVertexTotDof = 0, numCellTotDof = 0, numFaceTotDof[3] = {0, 0, 0}; 231 PetscSF sf; 232 PetscMPIInt rank; 233 const PetscMPIInt *neighbors; 234 PetscInt *localPoints; 235 PetscSFNode *remotePoints; 236 PetscInt nleaves = 0, nleavesCheck = 0; 237 PetscInt nC, nVx, nVy, nVz, nV, nxF, nXF, nyF, nYF, nzF, nZF; 238 PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart, zfEnd; 239 PetscInt f, v, c, xf, yf, zf, xn, yn, zn; 240 PetscErrorCode ierr; 241 242 PetscFunctionBegin; 243 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 244 ierr = MPI_Comm_rank(((PetscObject) dm)->comm, &rank);CHKERRQ(ierr); 245 ierr = DMDAGetNumCells(dm, &nC);CHKERRQ(ierr); 246 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr); 247 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 248 ierr = DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);CHKERRQ(ierr); 249 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 250 ierr = DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 251 ierr = DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);CHKERRQ(ierr); 252 xfStart = vEnd; xfEnd = xfStart+nXF; 253 yfStart = xfEnd; yfEnd = yfStart+nYF; 254 zfStart = yfEnd; zfEnd = zfStart+nZF; 255 if (zfEnd != fEnd) SETERRQ2(((PetscObject) dm)->comm, PETSC_ERR_PLIB, "Invalid face end %d, should be %d", zfEnd, fEnd); 256 /* Create local section */ 257 ierr = DMDAGetInfo(dm, 0,0,0,0,0,0,0, &numFields, 0,0,0,0,0);CHKERRQ(ierr); 258 for(f = 0; f < numFields; ++f) { 259 if (numVertexDof) {numVertexTotDof += numVertexDof[f];} 260 if (numCellDof) {numCellTotDof += numCellDof[f];} 261 if (numFaceDof) {numFaceTotDof[0] += numFaceDof[f*dim+0]; 262 numFaceTotDof[1] += dim > 1 ? numFaceDof[f*dim+1] : 0; 263 numFaceTotDof[2] += dim > 2 ? numFaceDof[f*dim+2] : 0;} 264 } 265 ierr = PetscSectionCreate(((PetscObject) dm)->comm, &dm->defaultSection);CHKERRQ(ierr); 266 if (numFields > 1) { 267 ierr = PetscSectionSetNumFields(dm->defaultSection, numFields);CHKERRQ(ierr); 268 for(f = 0; f < numFields; ++f) { 269 const char *name; 270 271 ierr = DMDAGetFieldName(dm, f, &name);CHKERRQ(ierr); 272 ierr = PetscSectionSetFieldName(dm->defaultSection, f, name);CHKERRQ(ierr); 273 if (numComp) { 274 ierr = PetscSectionSetFieldComponents(dm->defaultSection, f, numComp[f]);CHKERRQ(ierr); 275 } 276 } 277 } else { 278 numFields = 0; 279 } 280 ierr = PetscSectionSetChart(dm->defaultSection, pStart, pEnd);CHKERRQ(ierr); 281 if (numVertexDof) { 282 for(v = vStart; v < vEnd; ++v) { 283 for(f = 0; f < numFields; ++f) { 284 ierr = PetscSectionSetFieldDof(dm->defaultSection, v, f, numVertexDof[f]);CHKERRQ(ierr); 285 } 286 ierr = PetscSectionSetDof(dm->defaultSection, v, numVertexTotDof);CHKERRQ(ierr); 287 } 288 } 289 if (numFaceDof) { 290 for(xf = xfStart; xf < xfEnd; ++xf) { 291 for(f = 0; f < numFields; ++f) { 292 ierr = PetscSectionSetFieldDof(dm->defaultSection, xf, f, numFaceDof[f*dim+0]);CHKERRQ(ierr); 293 } 294 ierr = PetscSectionSetDof(dm->defaultSection, xf, numFaceTotDof[0]);CHKERRQ(ierr); 295 } 296 for(yf = yfStart; yf < yfEnd; ++yf) { 297 for(f = 0; f < numFields; ++f) { 298 ierr = PetscSectionSetFieldDof(dm->defaultSection, yf, f, numFaceDof[f*dim+1]);CHKERRQ(ierr); 299 } 300 ierr = PetscSectionSetDof(dm->defaultSection, yf, numFaceTotDof[1]);CHKERRQ(ierr); 301 } 302 for(zf = zfStart; zf < zfEnd; ++zf) { 303 for(f = 0; f < numFields; ++f) { 304 ierr = PetscSectionSetFieldDof(dm->defaultSection, zf, f, numFaceDof[f*dim+2]);CHKERRQ(ierr); 305 } 306 ierr = PetscSectionSetDof(dm->defaultSection, zf, numFaceTotDof[2]);CHKERRQ(ierr); 307 } 308 } 309 if (numCellDof) { 310 for(c = cStart; c < cEnd; ++c) { 311 for(f = 0; f < numFields; ++f) { 312 ierr = PetscSectionSetFieldDof(dm->defaultSection, c, f, numCellDof[f]);CHKERRQ(ierr); 313 } 314 ierr = PetscSectionSetDof(dm->defaultSection, c, numCellTotDof);CHKERRQ(ierr); 315 } 316 } 317 ierr = PetscSectionSetUp(dm->defaultSection);CHKERRQ(ierr); 318 /* Create mesh point SF */ 319 ierr = DMDAGetNeighbors(dm, &neighbors);CHKERRQ(ierr); 320 for(zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) { 321 for(yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) { 322 for(xn = 0; xn < 3; ++xn) { 323 const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0; 324 const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn]; 325 326 if (neighbor >= 0 && neighbor != rank) { 327 nleaves += (!xp ? nVx : 1) * (!yp ? nVy : 1) * (!zp ? nVz : 1); /* vertices */ 328 if (xp && !yp && !zp) { 329 nleaves += nxF; /* x faces */ 330 } else if (yp && !zp && !xp) { 331 nleaves += nyF; /* y faces */ 332 } else if (zp && !xp && !yp) { 333 nleaves += nzF; /* z faces */ 334 } 335 } 336 } 337 } 338 } 339 ierr = PetscMalloc2(nleaves,PetscInt,&localPoints,nleaves,PetscSFNode,&remotePoints);CHKERRQ(ierr); 340 for(zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) { 341 for(yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) { 342 for(xn = 0; xn < 3; ++xn) { 343 const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0; 344 const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn]; 345 346 if (neighbor >= 0 && neighbor != rank) { 347 if (xp < 0) { /* left */ 348 if (yp < 0) { /* bottom */ 349 if (zp < 0) { /* back */ 350 nleavesCheck += 1; /* left bottom back vertex */ 351 } else if (zp > 0) { /* front */ 352 nleavesCheck += 1; /* left bottom front vertex */ 353 } else { 354 nleavesCheck += nVz; /* left bottom vertices */ 355 } 356 } else if (yp > 0) { /* top */ 357 if (zp < 0) { /* back */ 358 nleavesCheck += 1; /* left top back vertex */ 359 } else if (zp > 0) { /* front */ 360 nleavesCheck += 1; /* left top front vertex */ 361 } else { 362 nleavesCheck += nVz; /* left top vertices */ 363 } 364 } else { 365 if (zp < 0) { /* back */ 366 nleavesCheck += nVy; /* left back vertices */ 367 } else if (zp > 0) { /* front */ 368 nleavesCheck += nVy; /* left front vertices */ 369 } else { 370 nleavesCheck += nVy*nVz; /* left vertices */ 371 nleavesCheck += nxF; /* left faces */ 372 } 373 } 374 } else if (xp > 0) { /* right */ 375 if (yp < 0) { /* bottom */ 376 if (zp < 0) { /* back */ 377 nleavesCheck += 1; /* right bottom back vertex */ 378 } else if (zp > 0) { /* front */ 379 nleavesCheck += 1; /* right bottom front vertex */ 380 } else { 381 nleavesCheck += nVz; /* right bottom vertices */ 382 } 383 } else if (yp > 0) { /* top */ 384 if (zp < 0) { /* back */ 385 nleavesCheck += 1; /* right top back vertex */ 386 } else if (zp > 0) { /* front */ 387 nleavesCheck += 1; /* right top front vertex */ 388 } else { 389 nleavesCheck += nVz; /* right top vertices */ 390 } 391 } else { 392 if (zp < 0) { /* back */ 393 nleavesCheck += nVy; /* right back vertices */ 394 } else if (zp > 0) { /* front */ 395 nleavesCheck += nVy; /* right front vertices */ 396 } else { 397 nleavesCheck += nVy*nVz; /* right vertices */ 398 nleavesCheck += nxF; /* right faces */ 399 } 400 } 401 } else { 402 if (yp < 0) { /* bottom */ 403 if (zp < 0) { /* back */ 404 nleavesCheck += nVx; /* bottom back vertices */ 405 } else if (zp > 0) { /* front */ 406 nleavesCheck += nVx; /* bottom front vertices */ 407 } else { 408 nleavesCheck += nVx*nVz; /* bottom vertices */ 409 nleavesCheck += nyF; /* bottom faces */ 410 } 411 } else if (yp > 0) { /* top */ 412 if (zp < 0) { /* back */ 413 nleavesCheck += nVx; /* top back vertices */ 414 } else if (zp > 0) { /* front */ 415 nleavesCheck += nVx; /* top front vertices */ 416 } else { 417 nleavesCheck += nVx*nVz; /* top vertices */ 418 nleavesCheck += nyF; /* top faces */ 419 } 420 } else { 421 if (zp < 0) { /* back */ 422 nleavesCheck += nVx*nVy; /* back vertices */ 423 nleavesCheck += nzF; /* back faces */ 424 } else if (zp > 0) { /* front */ 425 nleavesCheck += nVx*nVy; /* front vertices */ 426 nleavesCheck += nzF; /* front faces */ 427 } else { 428 /* Nothing is shared from the interior */ 429 } 430 } 431 } 432 } 433 } 434 } 435 } 436 if (nleaves != nleavesCheck) SETERRQ2(((PetscObject) dm)->comm, PETSC_ERR_PLIB, "The number of leaves %d did not match the number of remote leaves %d", nleaves, nleavesCheck); 437 ierr = PetscSFCreate(((PetscObject) dm)->comm, &sf);CHKERRQ(ierr); 438 ierr = PetscSFSetGraph(sf, pEnd, nleaves, localPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 439 /* Create global section */ 440 ierr = PetscSectionCreateGlobalSection(dm->defaultSection, sf, &dm->defaultGlobalSection);CHKERRQ(ierr); 441 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 442 /* Create default SF */ 443 ierr = DMCreateDefaultSF(dm, dm->defaultSection, dm->defaultGlobalSection);CHKERRQ(ierr); 444 PetscFunctionReturn(0); 445 } 446 447 /* ------------------------------------------------------------------- */ 448 #if defined(PETSC_HAVE_ADIC) 449 450 EXTERN_C_BEGIN 451 #include <adic/ad_utils.h> 452 EXTERN_C_END 453 454 #undef __FUNCT__ 455 #define __FUNCT__ "DMDAGetAdicArray" 456 /*@C 457 DMDAGetAdicArray - Gets an array of derivative types for a DMDA 458 459 Input Parameter: 460 + da - information about my local patch 461 - ghosted - do you want arrays for the ghosted or nonghosted patch 462 463 Output Parameters: 464 + vptr - array data structured to be passed to ad_FormFunctionLocal() 465 . array_start - actual start of 1d array of all values that adiC can access directly (may be null) 466 - tdof - total number of degrees of freedom represented in array_start (may be null) 467 468 Notes: 469 The vector values are NOT initialized and may have garbage in them, so you may need 470 to zero them. 471 472 Returns the same type of object as the DMDAVecGetArray() except its elements are 473 derivative types instead of PetscScalars 474 475 Level: advanced 476 477 .seealso: DMDARestoreAdicArray() 478 479 @*/ 480 PetscErrorCode DMDAGetAdicArray(DM da,PetscBool ghosted,void *vptr,void *array_start,PetscInt *tdof) 481 { 482 PetscErrorCode ierr; 483 PetscInt j,i,deriv_type_size,xs,ys,xm,ym,zs,zm,itdof; 484 char *iarray_start; 485 void **iptr = (void**)vptr; 486 DM_DA *dd = (DM_DA*)da->data; 487 488 PetscFunctionBegin; 489 PetscValidHeaderSpecific(da,DM_CLASSID,1); 490 if (ghosted) { 491 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 492 if (dd->adarrayghostedin[i]) { 493 *iptr = dd->adarrayghostedin[i]; 494 iarray_start = (char*)dd->adstartghostedin[i]; 495 itdof = dd->ghostedtdof; 496 dd->adarrayghostedin[i] = PETSC_NULL; 497 dd->adstartghostedin[i] = PETSC_NULL; 498 499 goto done; 500 } 501 } 502 xs = dd->Xs; 503 ys = dd->Ys; 504 zs = dd->Zs; 505 xm = dd->Xe-dd->Xs; 506 ym = dd->Ye-dd->Ys; 507 zm = dd->Ze-dd->Zs; 508 } else { 509 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 510 if (dd->adarrayin[i]) { 511 *iptr = dd->adarrayin[i]; 512 iarray_start = (char*)dd->adstartin[i]; 513 itdof = dd->tdof; 514 dd->adarrayin[i] = PETSC_NULL; 515 dd->adstartin[i] = PETSC_NULL; 516 517 goto done; 518 } 519 } 520 xs = dd->xs; 521 ys = dd->ys; 522 zs = dd->zs; 523 xm = dd->xe-dd->xs; 524 ym = dd->ye-dd->ys; 525 zm = dd->ze-dd->zs; 526 } 527 deriv_type_size = PetscADGetDerivTypeSize(); 528 529 switch (dd->dim) { 530 case 1: { 531 void *ptr; 532 itdof = xm; 533 534 ierr = PetscMalloc(xm*deriv_type_size,&iarray_start);CHKERRQ(ierr); 535 536 ptr = (void*)(iarray_start - xs*deriv_type_size); 537 *iptr = (void*)ptr; 538 break;} 539 case 2: { 540 void **ptr; 541 itdof = xm*ym; 542 543 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*deriv_type_size,&iarray_start);CHKERRQ(ierr); 544 545 ptr = (void**)(iarray_start + xm*ym*deriv_type_size - ys*sizeof(void*)); 546 for(j=ys;j<ys+ym;j++) { 547 ptr[j] = iarray_start + deriv_type_size*(xm*(j-ys) - xs); 548 } 549 *iptr = (void*)ptr; 550 break;} 551 case 3: { 552 void ***ptr,**bptr; 553 itdof = xm*ym*zm; 554 555 ierr = PetscMalloc((zm+1)*sizeof(void **)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*deriv_type_size,&iarray_start);CHKERRQ(ierr); 556 557 ptr = (void***)(iarray_start + xm*ym*zm*deriv_type_size - zs*sizeof(void*)); 558 bptr = (void**)(iarray_start + xm*ym*zm*deriv_type_size + zm*sizeof(void**)); 559 560 for(i=zs;i<zs+zm;i++) { 561 ptr[i] = bptr + ((i-zs)*ym - ys); 562 } 563 for (i=zs; i<zs+zm; i++) { 564 for (j=ys; j<ys+ym; j++) { 565 ptr[i][j] = iarray_start + deriv_type_size*(xm*ym*(i-zs) + xm*(j-ys) - xs); 566 } 567 } 568 569 *iptr = (void*)ptr; 570 break;} 571 default: 572 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 573 } 574 575 done: 576 /* add arrays to the checked out list */ 577 if (ghosted) { 578 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 579 if (!dd->adarrayghostedout[i]) { 580 dd->adarrayghostedout[i] = *iptr ; 581 dd->adstartghostedout[i] = iarray_start; 582 dd->ghostedtdof = itdof; 583 break; 584 } 585 } 586 } else { 587 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 588 if (!dd->adarrayout[i]) { 589 dd->adarrayout[i] = *iptr ; 590 dd->adstartout[i] = iarray_start; 591 dd->tdof = itdof; 592 break; 593 } 594 } 595 } 596 if (i == DMDA_MAX_AD_ARRAYS+1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Too many DMDA ADIC arrays obtained"); 597 if (tdof) *tdof = itdof; 598 if (array_start) *(void**)array_start = iarray_start; 599 PetscFunctionReturn(0); 600 } 601 602 #undef __FUNCT__ 603 #define __FUNCT__ "DMDARestoreAdicArray" 604 /*@C 605 DMDARestoreAdicArray - Restores an array of derivative types for a DMDA 606 607 Input Parameter: 608 + da - information about my local patch 609 - ghosted - do you want arrays for the ghosted or nonghosted patch 610 611 Output Parameters: 612 + ptr - array data structured to be passed to ad_FormFunctionLocal() 613 . array_start - actual start of 1d array of all values that adiC can access directly 614 - tdof - total number of degrees of freedom represented in array_start 615 616 Level: advanced 617 618 .seealso: DMDAGetAdicArray() 619 620 @*/ 621 PetscErrorCode DMDARestoreAdicArray(DM da,PetscBool ghosted,void *ptr,void *array_start,PetscInt *tdof) 622 { 623 PetscInt i; 624 void **iptr = (void**)ptr,iarray_start = 0; 625 626 PetscFunctionBegin; 627 PetscValidHeaderSpecific(da,DM_CLASSID,1); 628 if (ghosted) { 629 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 630 if (dd->adarrayghostedout[i] == *iptr) { 631 iarray_start = dd->adstartghostedout[i]; 632 dd->adarrayghostedout[i] = PETSC_NULL; 633 dd->adstartghostedout[i] = PETSC_NULL; 634 break; 635 } 636 } 637 if (!iarray_start) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Could not find array in checkout list"); 638 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 639 if (!dd->adarrayghostedin[i]){ 640 dd->adarrayghostedin[i] = *iptr; 641 dd->adstartghostedin[i] = iarray_start; 642 break; 643 } 644 } 645 } else { 646 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 647 if (dd->adarrayout[i] == *iptr) { 648 iarray_start = dd->adstartout[i]; 649 dd->adarrayout[i] = PETSC_NULL; 650 dd->adstartout[i] = PETSC_NULL; 651 break; 652 } 653 } 654 if (!iarray_start) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Could not find array in checkout list"); 655 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 656 if (!dd->adarrayin[i]){ 657 dd->adarrayin[i] = *iptr; 658 dd->adstartin[i] = iarray_start; 659 break; 660 } 661 } 662 } 663 PetscFunctionReturn(0); 664 } 665 666 #undef __FUNCT__ 667 #define __FUNCT__ "ad_DAGetArray" 668 PetscErrorCode ad_DAGetArray(DM da,PetscBool ghosted,void *iptr) 669 { 670 PetscErrorCode ierr; 671 PetscFunctionBegin; 672 ierr = DMDAGetAdicArray(da,ghosted,iptr,0,0);CHKERRQ(ierr); 673 PetscFunctionReturn(0); 674 } 675 676 #undef __FUNCT__ 677 #define __FUNCT__ "ad_DARestoreArray" 678 PetscErrorCode ad_DARestoreArray(DM da,PetscBool ghosted,void *iptr) 679 { 680 PetscErrorCode ierr; 681 PetscFunctionBegin; 682 ierr = DMDARestoreAdicArray(da,ghosted,iptr,0,0);CHKERRQ(ierr); 683 PetscFunctionReturn(0); 684 } 685 686 #endif 687 688 #undef __FUNCT__ 689 #define __FUNCT__ "DMDAGetArray" 690 /*@C 691 DMDAGetArray - Gets a work array for a DMDA 692 693 Input Parameter: 694 + da - information about my local patch 695 - ghosted - do you want arrays for the ghosted or nonghosted patch 696 697 Output Parameters: 698 . vptr - array data structured 699 700 Note: The vector values are NOT initialized and may have garbage in them, so you may need 701 to zero them. 702 703 Level: advanced 704 705 .seealso: DMDARestoreArray(), DMDAGetAdicArray() 706 707 @*/ 708 PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr) 709 { 710 PetscErrorCode ierr; 711 PetscInt j,i,xs,ys,xm,ym,zs,zm; 712 char *iarray_start; 713 void **iptr = (void**)vptr; 714 DM_DA *dd = (DM_DA*)da->data; 715 716 PetscFunctionBegin; 717 PetscValidHeaderSpecific(da,DM_CLASSID,1); 718 if (ghosted) { 719 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 720 if (dd->arrayghostedin[i]) { 721 *iptr = dd->arrayghostedin[i]; 722 iarray_start = (char*)dd->startghostedin[i]; 723 dd->arrayghostedin[i] = PETSC_NULL; 724 dd->startghostedin[i] = PETSC_NULL; 725 726 goto done; 727 } 728 } 729 xs = dd->Xs; 730 ys = dd->Ys; 731 zs = dd->Zs; 732 xm = dd->Xe-dd->Xs; 733 ym = dd->Ye-dd->Ys; 734 zm = dd->Ze-dd->Zs; 735 } else { 736 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 737 if (dd->arrayin[i]) { 738 *iptr = dd->arrayin[i]; 739 iarray_start = (char*)dd->startin[i]; 740 dd->arrayin[i] = PETSC_NULL; 741 dd->startin[i] = PETSC_NULL; 742 743 goto done; 744 } 745 } 746 xs = dd->xs; 747 ys = dd->ys; 748 zs = dd->zs; 749 xm = dd->xe-dd->xs; 750 ym = dd->ye-dd->ys; 751 zm = dd->ze-dd->zs; 752 } 753 754 switch (dd->dim) { 755 case 1: { 756 void *ptr; 757 758 ierr = PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 759 760 ptr = (void*)(iarray_start - xs*sizeof(PetscScalar)); 761 *iptr = (void*)ptr; 762 break;} 763 case 2: { 764 void **ptr; 765 766 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 767 768 ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*)); 769 for(j=ys;j<ys+ym;j++) { 770 ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs); 771 } 772 *iptr = (void*)ptr; 773 break;} 774 case 3: { 775 void ***ptr,**bptr; 776 777 ierr = PetscMalloc((zm+1)*sizeof(void **)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 778 779 ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*)); 780 bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**)); 781 for(i=zs;i<zs+zm;i++) { 782 ptr[i] = bptr + ((i-zs)*ym - ys); 783 } 784 for (i=zs; i<zs+zm; i++) { 785 for (j=ys; j<ys+ym; j++) { 786 ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs); 787 } 788 } 789 790 *iptr = (void*)ptr; 791 break;} 792 default: 793 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 794 } 795 796 done: 797 /* add arrays to the checked out list */ 798 if (ghosted) { 799 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 800 if (!dd->arrayghostedout[i]) { 801 dd->arrayghostedout[i] = *iptr ; 802 dd->startghostedout[i] = iarray_start; 803 break; 804 } 805 } 806 } else { 807 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 808 if (!dd->arrayout[i]) { 809 dd->arrayout[i] = *iptr ; 810 dd->startout[i] = iarray_start; 811 break; 812 } 813 } 814 } 815 PetscFunctionReturn(0); 816 } 817 818 #undef __FUNCT__ 819 #define __FUNCT__ "DMDARestoreArray" 820 /*@C 821 DMDARestoreArray - Restores an array of derivative types for a DMDA 822 823 Input Parameter: 824 + da - information about my local patch 825 . ghosted - do you want arrays for the ghosted or nonghosted patch 826 - vptr - array data structured to be passed to ad_FormFunctionLocal() 827 828 Level: advanced 829 830 .seealso: DMDAGetArray(), DMDAGetAdicArray() 831 832 @*/ 833 PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr) 834 { 835 PetscInt i; 836 void **iptr = (void**)vptr,*iarray_start = 0; 837 DM_DA *dd = (DM_DA*)da->data; 838 839 PetscFunctionBegin; 840 PetscValidHeaderSpecific(da,DM_CLASSID,1); 841 if (ghosted) { 842 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 843 if (dd->arrayghostedout[i] == *iptr) { 844 iarray_start = dd->startghostedout[i]; 845 dd->arrayghostedout[i] = PETSC_NULL; 846 dd->startghostedout[i] = PETSC_NULL; 847 break; 848 } 849 } 850 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 851 if (!dd->arrayghostedin[i]){ 852 dd->arrayghostedin[i] = *iptr; 853 dd->startghostedin[i] = iarray_start; 854 break; 855 } 856 } 857 } else { 858 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 859 if (dd->arrayout[i] == *iptr) { 860 iarray_start = dd->startout[i]; 861 dd->arrayout[i] = PETSC_NULL; 862 dd->startout[i] = PETSC_NULL; 863 break; 864 } 865 } 866 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 867 if (!dd->arrayin[i]){ 868 dd->arrayin[i] = *iptr; 869 dd->startin[i] = iarray_start; 870 break; 871 } 872 } 873 } 874 PetscFunctionReturn(0); 875 } 876 877 #undef __FUNCT__ 878 #define __FUNCT__ "DMDAGetAdicMFArray" 879 /*@C 880 DMDAGetAdicMFArray - Gets an array of derivative types for a DMDA for matrix-free ADIC. 881 882 Input Parameter: 883 + da - information about my local patch 884 - ghosted - do you want arrays for the ghosted or nonghosted patch? 885 886 Output Parameters: 887 + vptr - array data structured to be passed to ad_FormFunctionLocal() 888 . array_start - actual start of 1d array of all values that adiC can access directly (may be null) 889 - tdof - total number of degrees of freedom represented in array_start (may be null) 890 891 Notes: 892 The vector values are NOT initialized and may have garbage in them, so you may need 893 to zero them. 894 895 This routine returns the same type of object as the DMDAVecGetArray(), except its 896 elements are derivative types instead of PetscScalars. 897 898 Level: advanced 899 900 .seealso: DMDARestoreAdicMFArray(), DMDAGetArray(), DMDAGetAdicArray() 901 902 @*/ 903 PetscErrorCode DMDAGetAdicMFArray(DM da,PetscBool ghosted,void *vptr,void *array_start,PetscInt *tdof) 904 { 905 PetscErrorCode ierr; 906 PetscInt j,i,xs,ys,xm,ym,zs,zm,itdof = 0; 907 char *iarray_start; 908 void **iptr = (void**)vptr; 909 DM_DA *dd = (DM_DA*)da->data; 910 911 PetscFunctionBegin; 912 PetscValidHeaderSpecific(da,DM_CLASSID,1); 913 if (ghosted) { 914 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 915 if (dd->admfarrayghostedin[i]) { 916 *iptr = dd->admfarrayghostedin[i]; 917 iarray_start = (char*)dd->admfstartghostedin[i]; 918 itdof = dd->ghostedtdof; 919 dd->admfarrayghostedin[i] = PETSC_NULL; 920 dd->admfstartghostedin[i] = PETSC_NULL; 921 922 goto done; 923 } 924 } 925 xs = dd->Xs; 926 ys = dd->Ys; 927 zs = dd->Zs; 928 xm = dd->Xe-dd->Xs; 929 ym = dd->Ye-dd->Ys; 930 zm = dd->Ze-dd->Zs; 931 } else { 932 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 933 if (dd->admfarrayin[i]) { 934 *iptr = dd->admfarrayin[i]; 935 iarray_start = (char*)dd->admfstartin[i]; 936 itdof = dd->tdof; 937 dd->admfarrayin[i] = PETSC_NULL; 938 dd->admfstartin[i] = PETSC_NULL; 939 940 goto done; 941 } 942 } 943 xs = dd->xs; 944 ys = dd->ys; 945 zs = dd->zs; 946 xm = dd->xe-dd->xs; 947 ym = dd->ye-dd->ys; 948 zm = dd->ze-dd->zs; 949 } 950 951 switch (dd->dim) { 952 case 1: { 953 void *ptr; 954 itdof = xm; 955 956 ierr = PetscMalloc(xm*2*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 957 958 ptr = (void*)(iarray_start - xs*2*sizeof(PetscScalar)); 959 *iptr = (void*)ptr; 960 break;} 961 case 2: { 962 void **ptr; 963 itdof = xm*ym; 964 965 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*2*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 966 967 ptr = (void**)(iarray_start + xm*ym*2*sizeof(PetscScalar) - ys*sizeof(void*)); 968 for(j=ys;j<ys+ym;j++) { 969 ptr[j] = iarray_start + 2*sizeof(PetscScalar)*(xm*(j-ys) - xs); 970 } 971 *iptr = (void*)ptr; 972 break;} 973 case 3: { 974 void ***ptr,**bptr; 975 itdof = xm*ym*zm; 976 977 ierr = PetscMalloc((zm+1)*sizeof(void **)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*2*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 978 979 ptr = (void***)(iarray_start + xm*ym*zm*2*sizeof(PetscScalar) - zs*sizeof(void*)); 980 bptr = (void**)(iarray_start + xm*ym*zm*2*sizeof(PetscScalar) + zm*sizeof(void**)); 981 for(i=zs;i<zs+zm;i++) { 982 ptr[i] = bptr + ((i-zs)*ym* - ys)*sizeof(void*); 983 } 984 for (i=zs; i<zs+zm; i++) { 985 for (j=ys; j<ys+ym; j++) { 986 ptr[i][j] = iarray_start + 2*sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs); 987 } 988 } 989 990 *iptr = (void*)ptr; 991 break;} 992 default: 993 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 994 } 995 996 done: 997 /* add arrays to the checked out list */ 998 if (ghosted) { 999 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1000 if (!dd->admfarrayghostedout[i]) { 1001 dd->admfarrayghostedout[i] = *iptr ; 1002 dd->admfstartghostedout[i] = iarray_start; 1003 dd->ghostedtdof = itdof; 1004 break; 1005 } 1006 } 1007 } else { 1008 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1009 if (!dd->admfarrayout[i]) { 1010 dd->admfarrayout[i] = *iptr ; 1011 dd->admfstartout[i] = iarray_start; 1012 dd->tdof = itdof; 1013 break; 1014 } 1015 } 1016 } 1017 if (i == DMDA_MAX_AD_ARRAYS+1) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Too many DMDA ADIC arrays obtained"); 1018 if (tdof) *tdof = itdof; 1019 if (array_start) *(void**)array_start = iarray_start; 1020 PetscFunctionReturn(0); 1021 } 1022 1023 #undef __FUNCT__ 1024 #define __FUNCT__ "DMDAGetAdicMFArray4" 1025 PetscErrorCode DMDAGetAdicMFArray4(DM da,PetscBool ghosted,void *vptr,void *array_start,PetscInt *tdof) 1026 { 1027 PetscErrorCode ierr; 1028 PetscInt j,i,xs,ys,xm,ym,itdof = 0; 1029 char *iarray_start; 1030 void **iptr = (void**)vptr; 1031 DM_DA *dd = (DM_DA*)da->data; 1032 1033 PetscFunctionBegin; 1034 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1035 if (ghosted) { 1036 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1037 if (dd->admfarrayghostedin[i]) { 1038 *iptr = dd->admfarrayghostedin[i]; 1039 iarray_start = (char*)dd->admfstartghostedin[i]; 1040 itdof = dd->ghostedtdof; 1041 dd->admfarrayghostedin[i] = PETSC_NULL; 1042 dd->admfstartghostedin[i] = PETSC_NULL; 1043 1044 goto done; 1045 } 1046 } 1047 xs = dd->Xs; 1048 ys = dd->Ys; 1049 xm = dd->Xe-dd->Xs; 1050 ym = dd->Ye-dd->Ys; 1051 } else { 1052 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1053 if (dd->admfarrayin[i]) { 1054 *iptr = dd->admfarrayin[i]; 1055 iarray_start = (char*)dd->admfstartin[i]; 1056 itdof = dd->tdof; 1057 dd->admfarrayin[i] = PETSC_NULL; 1058 dd->admfstartin[i] = PETSC_NULL; 1059 1060 goto done; 1061 } 1062 } 1063 xs = dd->xs; 1064 ys = dd->ys; 1065 xm = dd->xe-dd->xs; 1066 ym = dd->ye-dd->ys; 1067 } 1068 1069 switch (dd->dim) { 1070 case 2: { 1071 void **ptr; 1072 itdof = xm*ym; 1073 1074 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*5*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1075 1076 ptr = (void**)(iarray_start + xm*ym*5*sizeof(PetscScalar) - ys*sizeof(void*)); 1077 for(j=ys;j<ys+ym;j++) { 1078 ptr[j] = iarray_start + 5*sizeof(PetscScalar)*(xm*(j-ys) - xs); 1079 } 1080 *iptr = (void*)ptr; 1081 break;} 1082 default: 1083 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 1084 } 1085 1086 done: 1087 /* add arrays to the checked out list */ 1088 if (ghosted) { 1089 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1090 if (!dd->admfarrayghostedout[i]) { 1091 dd->admfarrayghostedout[i] = *iptr ; 1092 dd->admfstartghostedout[i] = iarray_start; 1093 dd->ghostedtdof = itdof; 1094 break; 1095 } 1096 } 1097 } else { 1098 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1099 if (!dd->admfarrayout[i]) { 1100 dd->admfarrayout[i] = *iptr ; 1101 dd->admfstartout[i] = iarray_start; 1102 dd->tdof = itdof; 1103 break; 1104 } 1105 } 1106 } 1107 if (i == DMDA_MAX_AD_ARRAYS+1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Too many DMDA ADIC arrays obtained"); 1108 if (tdof) *tdof = itdof; 1109 if (array_start) *(void**)array_start = iarray_start; 1110 PetscFunctionReturn(0); 1111 } 1112 1113 #undef __FUNCT__ 1114 #define __FUNCT__ "DMDAGetAdicMFArray9" 1115 PetscErrorCode DMDAGetAdicMFArray9(DM da,PetscBool ghosted,void *vptr,void *array_start,PetscInt *tdof) 1116 { 1117 PetscErrorCode ierr; 1118 PetscInt j,i,xs,ys,xm,ym,itdof = 0; 1119 char *iarray_start; 1120 void **iptr = (void**)vptr; 1121 DM_DA *dd = (DM_DA*)da->data; 1122 1123 PetscFunctionBegin; 1124 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1125 if (ghosted) { 1126 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1127 if (dd->admfarrayghostedin[i]) { 1128 *iptr = dd->admfarrayghostedin[i]; 1129 iarray_start = (char*)dd->admfstartghostedin[i]; 1130 itdof = dd->ghostedtdof; 1131 dd->admfarrayghostedin[i] = PETSC_NULL; 1132 dd->admfstartghostedin[i] = PETSC_NULL; 1133 1134 goto done; 1135 } 1136 } 1137 xs = dd->Xs; 1138 ys = dd->Ys; 1139 xm = dd->Xe-dd->Xs; 1140 ym = dd->Ye-dd->Ys; 1141 } else { 1142 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1143 if (dd->admfarrayin[i]) { 1144 *iptr = dd->admfarrayin[i]; 1145 iarray_start = (char*)dd->admfstartin[i]; 1146 itdof = dd->tdof; 1147 dd->admfarrayin[i] = PETSC_NULL; 1148 dd->admfstartin[i] = PETSC_NULL; 1149 1150 goto done; 1151 } 1152 } 1153 xs = dd->xs; 1154 ys = dd->ys; 1155 xm = dd->xe-dd->xs; 1156 ym = dd->ye-dd->ys; 1157 } 1158 1159 switch (dd->dim) { 1160 case 2: { 1161 void **ptr; 1162 itdof = xm*ym; 1163 1164 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*10*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1165 1166 ptr = (void**)(iarray_start + xm*ym*10*sizeof(PetscScalar) - ys*sizeof(void*)); 1167 for(j=ys;j<ys+ym;j++) { 1168 ptr[j] = iarray_start + 10*sizeof(PetscScalar)*(xm*(j-ys) - xs); 1169 } 1170 *iptr = (void*)ptr; 1171 break;} 1172 default: 1173 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 1174 } 1175 1176 done: 1177 /* add arrays to the checked out list */ 1178 if (ghosted) { 1179 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1180 if (!dd->admfarrayghostedout[i]) { 1181 dd->admfarrayghostedout[i] = *iptr ; 1182 dd->admfstartghostedout[i] = iarray_start; 1183 dd->ghostedtdof = itdof; 1184 break; 1185 } 1186 } 1187 } else { 1188 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1189 if (!dd->admfarrayout[i]) { 1190 dd->admfarrayout[i] = *iptr ; 1191 dd->admfstartout[i] = iarray_start; 1192 dd->tdof = itdof; 1193 break; 1194 } 1195 } 1196 } 1197 if (i == DMDA_MAX_AD_ARRAYS+1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Too many DMDA ADIC arrays obtained"); 1198 if (tdof) *tdof = itdof; 1199 if (array_start) *(void**)array_start = iarray_start; 1200 PetscFunctionReturn(0); 1201 } 1202 1203 #undef __FUNCT__ 1204 #define __FUNCT__ "DMDAGetAdicMFArrayb" 1205 /*@C 1206 DMDAGetAdicMFArrayb - Gets an array of derivative types for a DMDA for matrix-free ADIC. 1207 1208 Input Parameter: 1209 + da - information about my local patch 1210 - ghosted - do you want arrays for the ghosted or nonghosted patch? 1211 1212 Output Parameters: 1213 + vptr - array data structured to be passed to ad_FormFunctionLocal() 1214 . array_start - actual start of 1d array of all values that adiC can access directly (may be null) 1215 - tdof - total number of degrees of freedom represented in array_start (may be null) 1216 1217 Notes: 1218 The vector values are NOT initialized and may have garbage in them, so you may need 1219 to zero them. 1220 1221 This routine returns the same type of object as the DMDAVecGetArray(), except its 1222 elements are derivative types instead of PetscScalars. 1223 1224 Level: advanced 1225 1226 .seealso: DMDARestoreAdicMFArray(), DMDAGetArray(), DMDAGetAdicArray() 1227 1228 @*/ 1229 PetscErrorCode DMDAGetAdicMFArrayb(DM da,PetscBool ghosted,void *vptr,void *array_start,PetscInt *tdof) 1230 { 1231 PetscErrorCode ierr; 1232 PetscInt j,i,xs,ys,xm,ym,zs,zm,itdof = 0; 1233 char *iarray_start; 1234 void **iptr = (void**)vptr; 1235 DM_DA *dd = (DM_DA*)da->data; 1236 PetscInt bs = dd->w,bs1 = bs+1; 1237 1238 PetscFunctionBegin; 1239 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1240 if (ghosted) { 1241 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1242 if (dd->admfarrayghostedin[i]) { 1243 *iptr = dd->admfarrayghostedin[i]; 1244 iarray_start = (char*)dd->admfstartghostedin[i]; 1245 itdof = dd->ghostedtdof; 1246 dd->admfarrayghostedin[i] = PETSC_NULL; 1247 dd->admfstartghostedin[i] = PETSC_NULL; 1248 1249 goto done; 1250 } 1251 } 1252 xs = dd->Xs; 1253 ys = dd->Ys; 1254 zs = dd->Zs; 1255 xm = dd->Xe-dd->Xs; 1256 ym = dd->Ye-dd->Ys; 1257 zm = dd->Ze-dd->Zs; 1258 } else { 1259 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1260 if (dd->admfarrayin[i]) { 1261 *iptr = dd->admfarrayin[i]; 1262 iarray_start = (char*)dd->admfstartin[i]; 1263 itdof = dd->tdof; 1264 dd->admfarrayin[i] = PETSC_NULL; 1265 dd->admfstartin[i] = PETSC_NULL; 1266 1267 goto done; 1268 } 1269 } 1270 xs = dd->xs; 1271 ys = dd->ys; 1272 zs = dd->zs; 1273 xm = dd->xe-dd->xs; 1274 ym = dd->ye-dd->ys; 1275 zm = dd->ze-dd->zs; 1276 } 1277 1278 switch (dd->dim) { 1279 case 1: { 1280 void *ptr; 1281 itdof = xm; 1282 1283 ierr = PetscMalloc(xm*bs1*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1284 1285 ptr = (void*)(iarray_start - xs*bs1*sizeof(PetscScalar)); 1286 *iptr = (void*)ptr; 1287 break;} 1288 case 2: { 1289 void **ptr; 1290 itdof = xm*ym; 1291 1292 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*bs1*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1293 1294 ptr = (void**)(iarray_start + xm*ym*bs1*sizeof(PetscScalar) - ys*sizeof(void*)); 1295 for(j=ys;j<ys+ym;j++) { 1296 ptr[j] = iarray_start + bs1*sizeof(PetscScalar)*(xm*(j-ys) - xs); 1297 } 1298 *iptr = (void*)ptr; 1299 break;} 1300 case 3: { 1301 void ***ptr,**bptr; 1302 itdof = xm*ym*zm; 1303 1304 ierr = PetscMalloc((zm+1)*sizeof(void **)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*bs1*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1305 1306 ptr = (void***)(iarray_start + xm*ym*zm*2*sizeof(PetscScalar) - zs*sizeof(void*)); 1307 bptr = (void**)(iarray_start + xm*ym*zm*2*sizeof(PetscScalar) + zm*sizeof(void**)); 1308 for(i=zs;i<zs+zm;i++) { 1309 ptr[i] = bptr + ((i-zs)*ym* - ys)*sizeof(void*); 1310 } 1311 for (i=zs; i<zs+zm; i++) { 1312 for (j=ys; j<ys+ym; j++) { 1313 ptr[i][j] = iarray_start + bs1*sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs); 1314 } 1315 } 1316 1317 *iptr = (void*)ptr; 1318 break;} 1319 default: 1320 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 1321 } 1322 1323 done: 1324 /* add arrays to the checked out list */ 1325 if (ghosted) { 1326 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1327 if (!dd->admfarrayghostedout[i]) { 1328 dd->admfarrayghostedout[i] = *iptr ; 1329 dd->admfstartghostedout[i] = iarray_start; 1330 dd->ghostedtdof = itdof; 1331 break; 1332 } 1333 } 1334 } else { 1335 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1336 if (!dd->admfarrayout[i]) { 1337 dd->admfarrayout[i] = *iptr ; 1338 dd->admfstartout[i] = iarray_start; 1339 dd->tdof = itdof; 1340 break; 1341 } 1342 } 1343 } 1344 if (i == DMDA_MAX_AD_ARRAYS+1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Too many DMDA ADIC arrays obtained"); 1345 if (tdof) *tdof = itdof; 1346 if (array_start) *(void**)array_start = iarray_start; 1347 PetscFunctionReturn(0); 1348 } 1349 1350 #undef __FUNCT__ 1351 #define __FUNCT__ "DMDARestoreAdicMFArray" 1352 /*@C 1353 DMDARestoreAdicMFArray - Restores an array of derivative types for a DMDA. 1354 1355 Input Parameter: 1356 + da - information about my local patch 1357 - ghosted - do you want arrays for the ghosted or nonghosted patch? 1358 1359 Output Parameters: 1360 + ptr - array data structure to be passed to ad_FormFunctionLocal() 1361 . array_start - actual start of 1d array of all values that adiC can access directly 1362 - tdof - total number of degrees of freedom represented in array_start 1363 1364 Level: advanced 1365 1366 .seealso: DMDAGetAdicArray() 1367 1368 @*/ 1369 PetscErrorCode DMDARestoreAdicMFArray(DM da,PetscBool ghosted,void *vptr,void *array_start,PetscInt *tdof) 1370 { 1371 PetscInt i; 1372 void **iptr = (void**)vptr,*iarray_start = 0; 1373 DM_DA *dd = (DM_DA*)da->data; 1374 1375 PetscFunctionBegin; 1376 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1377 if (ghosted) { 1378 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1379 if (dd->admfarrayghostedout[i] == *iptr) { 1380 iarray_start = dd->admfstartghostedout[i]; 1381 dd->admfarrayghostedout[i] = PETSC_NULL; 1382 dd->admfstartghostedout[i] = PETSC_NULL; 1383 break; 1384 } 1385 } 1386 if (!iarray_start) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Could not find array in checkout list"); 1387 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1388 if (!dd->admfarrayghostedin[i]){ 1389 dd->admfarrayghostedin[i] = *iptr; 1390 dd->admfstartghostedin[i] = iarray_start; 1391 break; 1392 } 1393 } 1394 } else { 1395 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1396 if (dd->admfarrayout[i] == *iptr) { 1397 iarray_start = dd->admfstartout[i]; 1398 dd->admfarrayout[i] = PETSC_NULL; 1399 dd->admfstartout[i] = PETSC_NULL; 1400 break; 1401 } 1402 } 1403 if (!iarray_start) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Could not find array in checkout list"); 1404 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1405 if (!dd->admfarrayin[i]){ 1406 dd->admfarrayin[i] = *iptr; 1407 dd->admfstartin[i] = iarray_start; 1408 break; 1409 } 1410 } 1411 } 1412 PetscFunctionReturn(0); 1413 } 1414 1415 #undef __FUNCT__ 1416 #define __FUNCT__ "admf_DAGetArray" 1417 PetscErrorCode admf_DAGetArray(DM da,PetscBool ghosted,void *iptr) 1418 { 1419 PetscErrorCode ierr; 1420 PetscFunctionBegin; 1421 ierr = DMDAGetAdicMFArray(da,ghosted,iptr,0,0);CHKERRQ(ierr); 1422 PetscFunctionReturn(0); 1423 } 1424 1425 #undef __FUNCT__ 1426 #define __FUNCT__ "admf_DARestoreArray" 1427 PetscErrorCode admf_DARestoreArray(DM da,PetscBool ghosted,void *iptr) 1428 { 1429 PetscErrorCode ierr; 1430 PetscFunctionBegin; 1431 ierr = DMDARestoreAdicMFArray(da,ghosted,iptr,0,0);CHKERRQ(ierr); 1432 PetscFunctionReturn(0); 1433 } 1434 1435