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