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