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, nL = 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 PetscInt xv, yv, zv; 346 347 if (neighbor >= 0 && neighbor != rank) { 348 if (xp < 0) { /* left */ 349 if (yp < 0) { /* bottom */ 350 if (zp < 0) { /* back */ 351 nleavesCheck += 1; /* left bottom back vertex */ 352 const PetscInt localVertex = ( 0*nVy + 0)*nVx + 0 + nC; 353 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 354 355 localPoints[nL] = localVertex; 356 remotePoints[nL].rank = neighbor; 357 remotePoints[nL].index = remoteVertex; 358 ++nL; 359 } else if (zp > 0) { /* front */ 360 nleavesCheck += 1; /* left bottom front vertex */ 361 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; 362 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 363 364 localPoints[nL] = localVertex; 365 remotePoints[nL].rank = neighbor; 366 remotePoints[nL].index = remoteVertex; 367 ++nL; 368 } else { 369 nleavesCheck += nVz; /* left bottom vertices */ 370 for(zv = 0; zv < nVz; ++zv, ++nL) { 371 const PetscInt localVertex = (zv*nVy + 0)*nVx + 0 + nC; 372 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 373 374 localPoints[nL] = localVertex; 375 remotePoints[nL].rank = neighbor; 376 remotePoints[nL].index = remoteVertex; 377 } 378 } 379 } else if (yp > 0) { /* top */ 380 if (zp < 0) { /* back */ 381 nleavesCheck += 1; /* left top back vertex */ 382 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; 383 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 384 385 localPoints[nL] = localVertex; 386 remotePoints[nL].rank = neighbor; 387 remotePoints[nL].index = remoteVertex; 388 ++nL; 389 } else if (zp > 0) { /* front */ 390 nleavesCheck += 1; /* left top front vertex */ 391 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; 392 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 393 394 localPoints[nL] = localVertex; 395 remotePoints[nL].rank = neighbor; 396 remotePoints[nL].index = remoteVertex; 397 ++nL; 398 } else { 399 nleavesCheck += nVz; /* left top vertices */ 400 for(zv = 0; zv < nVz; ++zv, ++nL) { 401 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; 402 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 403 404 localPoints[nL] = localVertex; 405 remotePoints[nL].rank = neighbor; 406 remotePoints[nL].index = remoteVertex; 407 } 408 } 409 } else { 410 if (zp < 0) { /* back */ 411 nleavesCheck += nVy; /* left back vertices */ 412 for(yv = 0; yv < nVy; ++yv, ++nL) { 413 const PetscInt localVertex = ( 0*nVy + yv)*nVx + 0 + nC; 414 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 415 416 localPoints[nL] = localVertex; 417 remotePoints[nL].rank = neighbor; 418 remotePoints[nL].index = remoteVertex; 419 } 420 } else if (zp > 0) { /* front */ 421 nleavesCheck += nVy; /* left front vertices */ 422 for(yv = 0; yv < nVy; ++yv, ++nL) { 423 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; 424 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 425 426 localPoints[nL] = localVertex; 427 remotePoints[nL].rank = neighbor; 428 remotePoints[nL].index = remoteVertex; 429 } 430 } else { 431 nleavesCheck += nVy*nVz; /* left vertices */ 432 for(zv = 0; zv < nVz; ++zv) { 433 for(yv = 0; yv < nVy; ++yv, ++nL) { 434 const PetscInt localVertex = (zv*nVy + yv)*nVx + 0 + nC; 435 const PetscInt remoteVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 436 437 localPoints[nL] = localVertex; 438 remotePoints[nL].rank = neighbor; 439 remotePoints[nL].index = remoteVertex; 440 } 441 } 442 nleavesCheck += nxF; /* left faces */ 443 for(xf = 0; xf < nxF; ++xf, ++nL) { 444 /* THIS IS WRONG */ 445 const PetscInt localFace = 0 + nC+nV; 446 const PetscInt remoteFace = 0 + nC+nV; 447 448 localPoints[nL] = localFace; 449 remotePoints[nL].rank = neighbor; 450 remotePoints[nL].index = remoteFace; 451 } 452 } 453 } 454 } else if (xp > 0) { /* right */ 455 if (yp < 0) { /* bottom */ 456 if (zp < 0) { /* back */ 457 nleavesCheck += 1; /* right bottom back vertex */ 458 const PetscInt localVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; 459 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 460 461 localPoints[nL] = localVertex; 462 remotePoints[nL].rank = neighbor; 463 remotePoints[nL].index = remoteVertex; 464 ++nL; 465 } else if (zp > 0) { /* front */ 466 nleavesCheck += 1; /* right bottom front vertex */ 467 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; 468 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 469 470 localPoints[nL] = localVertex; 471 remotePoints[nL].rank = neighbor; 472 remotePoints[nL].index = remoteVertex; 473 ++nL; 474 } else { 475 nleavesCheck += nVz; /* right bottom vertices */ 476 for(zv = 0; zv < nVz; ++zv, ++nL) { 477 const PetscInt localVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; 478 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 479 480 localPoints[nL] = localVertex; 481 remotePoints[nL].rank = neighbor; 482 remotePoints[nL].index = remoteVertex; 483 } 484 } 485 } else if (yp > 0) { /* top */ 486 if (zp < 0) { /* back */ 487 nleavesCheck += 1; /* right top back vertex */ 488 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; 489 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 490 491 localPoints[nL] = localVertex; 492 remotePoints[nL].rank = neighbor; 493 remotePoints[nL].index = remoteVertex; 494 ++nL; 495 } else if (zp > 0) { /* front */ 496 nleavesCheck += 1; /* right top front vertex */ 497 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; 498 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 499 500 localPoints[nL] = localVertex; 501 remotePoints[nL].rank = neighbor; 502 remotePoints[nL].index = remoteVertex; 503 ++nL; 504 } else { 505 nleavesCheck += nVz; /* right top vertices */ 506 for(zv = 0; zv < nVz; ++zv, ++nL) { 507 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; 508 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 509 510 localPoints[nL] = localVertex; 511 remotePoints[nL].rank = neighbor; 512 remotePoints[nL].index = remoteVertex; 513 } 514 } 515 } else { 516 if (zp < 0) { /* back */ 517 nleavesCheck += nVy; /* right back vertices */ 518 for(yv = 0; yv < nVy; ++yv, ++nL) { 519 const PetscInt localVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; 520 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 521 522 localPoints[nL] = localVertex; 523 remotePoints[nL].rank = neighbor; 524 remotePoints[nL].index = remoteVertex; 525 } 526 } else if (zp > 0) { /* front */ 527 nleavesCheck += nVy; /* right front vertices */ 528 for(yv = 0; yv < nVy; ++yv, ++nL) { 529 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; 530 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 531 532 localPoints[nL] = localVertex; 533 remotePoints[nL].rank = neighbor; 534 remotePoints[nL].index = remoteVertex; 535 } 536 } else { 537 nleavesCheck += nVy*nVz; /* right vertices */ 538 for(zv = 0; zv < nVz; ++zv) { 539 for(yv = 0; yv < nVy; ++yv, ++nL) { 540 const PetscInt localVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; 541 const PetscInt remoteVertex = (zv*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 542 543 localPoints[nL] = localVertex; 544 remotePoints[nL].rank = neighbor; 545 remotePoints[nL].index = remoteVertex; 546 } 547 } 548 nleavesCheck += nxF; /* right faces */ 549 for(xf = 0; xf < nxF; ++xf, ++nL) { 550 /* THIS IS WRONG */ 551 const PetscInt localFace = 0 + nC+nV; 552 const PetscInt remoteFace = 0 + nC+nV; 553 554 localPoints[nL] = localFace; 555 remotePoints[nL].rank = neighbor; 556 remotePoints[nL].index = remoteFace; 557 } 558 } 559 } 560 } else { 561 if (yp < 0) { /* bottom */ 562 if (zp < 0) { /* back */ 563 nleavesCheck += nVx; /* bottom back vertices */ 564 for(xv = 0; xv < nVx; ++xv, ++nL) { 565 const PetscInt localVertex = ( 0*nVy + 0)*nVx + xv + nC; 566 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 567 568 localPoints[nL] = localVertex; 569 remotePoints[nL].rank = neighbor; 570 remotePoints[nL].index = remoteVertex; 571 } 572 } else if (zp > 0) { /* front */ 573 nleavesCheck += nVx; /* bottom front vertices */ 574 for(xv = 0; xv < nVx; ++xv, ++nL) { 575 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; 576 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 577 578 localPoints[nL] = localVertex; 579 remotePoints[nL].rank = neighbor; 580 remotePoints[nL].index = remoteVertex; 581 } 582 } else { 583 nleavesCheck += nVx*nVz; /* bottom vertices */ 584 for(zv = 0; zv < nVz; ++zv) { 585 for(xv = 0; xv < nVx; ++xv, ++nL) { 586 const PetscInt localVertex = (zv*nVy + 0)*nVx + xv + nC; 587 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 588 589 localPoints[nL] = localVertex; 590 remotePoints[nL].rank = neighbor; 591 remotePoints[nL].index = remoteVertex; 592 } 593 } 594 nleavesCheck += nyF; /* bottom faces */ 595 for(yf = 0; yf < nyF; ++yf, ++nL) { 596 /* THIS IS WRONG */ 597 const PetscInt localFace = 0 + nC+nV; 598 const PetscInt remoteFace = 0 + nC+nV; 599 600 localPoints[nL] = localFace; 601 remotePoints[nL].rank = neighbor; 602 remotePoints[nL].index = remoteFace; 603 } 604 } 605 } else if (yp > 0) { /* top */ 606 if (zp < 0) { /* back */ 607 nleavesCheck += nVx; /* top back vertices */ 608 for(xv = 0; xv < nVx; ++xv, ++nL) { 609 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; 610 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 611 612 localPoints[nL] = localVertex; 613 remotePoints[nL].rank = neighbor; 614 remotePoints[nL].index = remoteVertex; 615 } 616 } else if (zp > 0) { /* front */ 617 nleavesCheck += nVx; /* top front vertices */ 618 for(xv = 0; xv < nVx; ++xv, ++nL) { 619 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; 620 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 621 622 localPoints[nL] = localVertex; 623 remotePoints[nL].rank = neighbor; 624 remotePoints[nL].index = remoteVertex; 625 } 626 } else { 627 nleavesCheck += nVx*nVz; /* top vertices */ 628 for(zv = 0; zv < nVz; ++zv) { 629 for(xv = 0; xv < nVx; ++xv, ++nL) { 630 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + xv + nC; 631 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 632 633 localPoints[nL] = localVertex; 634 remotePoints[nL].rank = neighbor; 635 remotePoints[nL].index = remoteVertex; 636 } 637 } 638 nleavesCheck += nyF; /* top faces */ 639 for(yf = 0; yf < nyF; ++yf, ++nL) { 640 /* THIS IS WRONG */ 641 const PetscInt localFace = 0 + nC+nV; 642 const PetscInt remoteFace = 0 + nC+nV; 643 644 localPoints[nL] = localFace; 645 remotePoints[nL].rank = neighbor; 646 remotePoints[nL].index = remoteFace; 647 } 648 } 649 } else { 650 if (zp < 0) { /* back */ 651 nleavesCheck += nVx*nVy; /* back vertices */ 652 for(yv = 0; yv < nVy; ++yv) { 653 for(xv = 0; xv < nVx; ++xv, ++nL) { 654 const PetscInt localVertex = ( 0*nVy + yv)*nVx + xv + nC; 655 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 656 657 localPoints[nL] = localVertex; 658 remotePoints[nL].rank = neighbor; 659 remotePoints[nL].index = remoteVertex; 660 } 661 } 662 nleavesCheck += nzF; /* back faces */ 663 for(zf = 0; zf < nzF; ++zf, ++nL) { 664 /* THIS IS WRONG */ 665 const PetscInt localFace = 0 + nC+nV; 666 const PetscInt remoteFace = 0 + nC+nV; 667 668 localPoints[nL] = localFace; 669 remotePoints[nL].rank = neighbor; 670 remotePoints[nL].index = remoteFace; 671 } 672 } else if (zp > 0) { /* front */ 673 nleavesCheck += nVx*nVy; /* front vertices */ 674 for(yv = 0; yv < nVy; ++yv) { 675 for(xv = 0; xv < nVx; ++xv, ++nL) { 676 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; 677 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 678 679 localPoints[nL] = localVertex; 680 remotePoints[nL].rank = neighbor; 681 remotePoints[nL].index = remoteVertex; 682 } 683 } 684 nleavesCheck += nzF; /* front faces */ 685 for(zf = 0; zf < nzF; ++zf, ++nL) { 686 /* THIS IS WRONG */ 687 const PetscInt localFace = 0 + nC+nV; 688 const PetscInt remoteFace = 0 + nC+nV; 689 690 localPoints[nL] = localFace; 691 remotePoints[nL].rank = neighbor; 692 remotePoints[nL].index = remoteFace; 693 } 694 } else { 695 /* Nothing is shared from the interior */ 696 } 697 } 698 } 699 } 700 } 701 } 702 } 703 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); 704 ierr = PetscSFCreate(((PetscObject) dm)->comm, &sf);CHKERRQ(ierr); 705 ierr = PetscSFSetGraph(sf, pEnd, nleaves, localPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 706 /* Create global section */ 707 ierr = PetscSectionCreateGlobalSection(dm->defaultSection, sf, &dm->defaultGlobalSection);CHKERRQ(ierr); 708 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 709 /* Create default SF */ 710 ierr = DMCreateDefaultSF(dm, dm->defaultSection, dm->defaultGlobalSection);CHKERRQ(ierr); 711 PetscFunctionReturn(0); 712 } 713 714 /* ------------------------------------------------------------------- */ 715 #if defined(PETSC_HAVE_ADIC) 716 717 EXTERN_C_BEGIN 718 #include <adic/ad_utils.h> 719 EXTERN_C_END 720 721 #undef __FUNCT__ 722 #define __FUNCT__ "DMDAGetAdicArray" 723 /*@C 724 DMDAGetAdicArray - Gets an array of derivative types for a DMDA 725 726 Input Parameter: 727 + da - information about my local patch 728 - ghosted - do you want arrays for the ghosted or nonghosted patch 729 730 Output Parameters: 731 + vptr - array data structured to be passed to ad_FormFunctionLocal() 732 . array_start - actual start of 1d array of all values that adiC can access directly (may be null) 733 - tdof - total number of degrees of freedom represented in array_start (may be null) 734 735 Notes: 736 The vector values are NOT initialized and may have garbage in them, so you may need 737 to zero them. 738 739 Returns the same type of object as the DMDAVecGetArray() except its elements are 740 derivative types instead of PetscScalars 741 742 Level: advanced 743 744 .seealso: DMDARestoreAdicArray() 745 746 @*/ 747 PetscErrorCode DMDAGetAdicArray(DM da,PetscBool ghosted,void *vptr,void *array_start,PetscInt *tdof) 748 { 749 PetscErrorCode ierr; 750 PetscInt j,i,deriv_type_size,xs,ys,xm,ym,zs,zm,itdof; 751 char *iarray_start; 752 void **iptr = (void**)vptr; 753 DM_DA *dd = (DM_DA*)da->data; 754 755 PetscFunctionBegin; 756 PetscValidHeaderSpecific(da,DM_CLASSID,1); 757 if (ghosted) { 758 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 759 if (dd->adarrayghostedin[i]) { 760 *iptr = dd->adarrayghostedin[i]; 761 iarray_start = (char*)dd->adstartghostedin[i]; 762 itdof = dd->ghostedtdof; 763 dd->adarrayghostedin[i] = PETSC_NULL; 764 dd->adstartghostedin[i] = PETSC_NULL; 765 766 goto done; 767 } 768 } 769 xs = dd->Xs; 770 ys = dd->Ys; 771 zs = dd->Zs; 772 xm = dd->Xe-dd->Xs; 773 ym = dd->Ye-dd->Ys; 774 zm = dd->Ze-dd->Zs; 775 } else { 776 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 777 if (dd->adarrayin[i]) { 778 *iptr = dd->adarrayin[i]; 779 iarray_start = (char*)dd->adstartin[i]; 780 itdof = dd->tdof; 781 dd->adarrayin[i] = PETSC_NULL; 782 dd->adstartin[i] = PETSC_NULL; 783 784 goto done; 785 } 786 } 787 xs = dd->xs; 788 ys = dd->ys; 789 zs = dd->zs; 790 xm = dd->xe-dd->xs; 791 ym = dd->ye-dd->ys; 792 zm = dd->ze-dd->zs; 793 } 794 deriv_type_size = PetscADGetDerivTypeSize(); 795 796 switch (dd->dim) { 797 case 1: { 798 void *ptr; 799 itdof = xm; 800 801 ierr = PetscMalloc(xm*deriv_type_size,&iarray_start);CHKERRQ(ierr); 802 803 ptr = (void*)(iarray_start - xs*deriv_type_size); 804 *iptr = (void*)ptr; 805 break;} 806 case 2: { 807 void **ptr; 808 itdof = xm*ym; 809 810 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*deriv_type_size,&iarray_start);CHKERRQ(ierr); 811 812 ptr = (void**)(iarray_start + xm*ym*deriv_type_size - ys*sizeof(void*)); 813 for(j=ys;j<ys+ym;j++) { 814 ptr[j] = iarray_start + deriv_type_size*(xm*(j-ys) - xs); 815 } 816 *iptr = (void*)ptr; 817 break;} 818 case 3: { 819 void ***ptr,**bptr; 820 itdof = xm*ym*zm; 821 822 ierr = PetscMalloc((zm+1)*sizeof(void **)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*deriv_type_size,&iarray_start);CHKERRQ(ierr); 823 824 ptr = (void***)(iarray_start + xm*ym*zm*deriv_type_size - zs*sizeof(void*)); 825 bptr = (void**)(iarray_start + xm*ym*zm*deriv_type_size + zm*sizeof(void**)); 826 827 for(i=zs;i<zs+zm;i++) { 828 ptr[i] = bptr + ((i-zs)*ym - ys); 829 } 830 for (i=zs; i<zs+zm; i++) { 831 for (j=ys; j<ys+ym; j++) { 832 ptr[i][j] = iarray_start + deriv_type_size*(xm*ym*(i-zs) + xm*(j-ys) - xs); 833 } 834 } 835 836 *iptr = (void*)ptr; 837 break;} 838 default: 839 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 840 } 841 842 done: 843 /* add arrays to the checked out list */ 844 if (ghosted) { 845 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 846 if (!dd->adarrayghostedout[i]) { 847 dd->adarrayghostedout[i] = *iptr ; 848 dd->adstartghostedout[i] = iarray_start; 849 dd->ghostedtdof = itdof; 850 break; 851 } 852 } 853 } else { 854 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 855 if (!dd->adarrayout[i]) { 856 dd->adarrayout[i] = *iptr ; 857 dd->adstartout[i] = iarray_start; 858 dd->tdof = itdof; 859 break; 860 } 861 } 862 } 863 if (i == DMDA_MAX_AD_ARRAYS+1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Too many DMDA ADIC arrays obtained"); 864 if (tdof) *tdof = itdof; 865 if (array_start) *(void**)array_start = iarray_start; 866 PetscFunctionReturn(0); 867 } 868 869 #undef __FUNCT__ 870 #define __FUNCT__ "DMDARestoreAdicArray" 871 /*@C 872 DMDARestoreAdicArray - Restores an array of derivative types for a DMDA 873 874 Input Parameter: 875 + da - information about my local patch 876 - ghosted - do you want arrays for the ghosted or nonghosted patch 877 878 Output Parameters: 879 + ptr - array data structured to be passed to ad_FormFunctionLocal() 880 . array_start - actual start of 1d array of all values that adiC can access directly 881 - tdof - total number of degrees of freedom represented in array_start 882 883 Level: advanced 884 885 .seealso: DMDAGetAdicArray() 886 887 @*/ 888 PetscErrorCode DMDARestoreAdicArray(DM da,PetscBool ghosted,void *ptr,void *array_start,PetscInt *tdof) 889 { 890 PetscInt i; 891 void **iptr = (void**)ptr,iarray_start = 0; 892 893 PetscFunctionBegin; 894 PetscValidHeaderSpecific(da,DM_CLASSID,1); 895 if (ghosted) { 896 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 897 if (dd->adarrayghostedout[i] == *iptr) { 898 iarray_start = dd->adstartghostedout[i]; 899 dd->adarrayghostedout[i] = PETSC_NULL; 900 dd->adstartghostedout[i] = PETSC_NULL; 901 break; 902 } 903 } 904 if (!iarray_start) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Could not find array in checkout list"); 905 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 906 if (!dd->adarrayghostedin[i]){ 907 dd->adarrayghostedin[i] = *iptr; 908 dd->adstartghostedin[i] = iarray_start; 909 break; 910 } 911 } 912 } else { 913 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 914 if (dd->adarrayout[i] == *iptr) { 915 iarray_start = dd->adstartout[i]; 916 dd->adarrayout[i] = PETSC_NULL; 917 dd->adstartout[i] = PETSC_NULL; 918 break; 919 } 920 } 921 if (!iarray_start) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Could not find array in checkout list"); 922 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 923 if (!dd->adarrayin[i]){ 924 dd->adarrayin[i] = *iptr; 925 dd->adstartin[i] = iarray_start; 926 break; 927 } 928 } 929 } 930 PetscFunctionReturn(0); 931 } 932 933 #undef __FUNCT__ 934 #define __FUNCT__ "ad_DAGetArray" 935 PetscErrorCode ad_DAGetArray(DM da,PetscBool ghosted,void *iptr) 936 { 937 PetscErrorCode ierr; 938 PetscFunctionBegin; 939 ierr = DMDAGetAdicArray(da,ghosted,iptr,0,0);CHKERRQ(ierr); 940 PetscFunctionReturn(0); 941 } 942 943 #undef __FUNCT__ 944 #define __FUNCT__ "ad_DARestoreArray" 945 PetscErrorCode ad_DARestoreArray(DM da,PetscBool ghosted,void *iptr) 946 { 947 PetscErrorCode ierr; 948 PetscFunctionBegin; 949 ierr = DMDARestoreAdicArray(da,ghosted,iptr,0,0);CHKERRQ(ierr); 950 PetscFunctionReturn(0); 951 } 952 953 #endif 954 955 #undef __FUNCT__ 956 #define __FUNCT__ "DMDAGetArray" 957 /*@C 958 DMDAGetArray - Gets a work array for a DMDA 959 960 Input Parameter: 961 + da - information about my local patch 962 - ghosted - do you want arrays for the ghosted or nonghosted patch 963 964 Output Parameters: 965 . vptr - array data structured 966 967 Note: The vector values are NOT initialized and may have garbage in them, so you may need 968 to zero them. 969 970 Level: advanced 971 972 .seealso: DMDARestoreArray(), DMDAGetAdicArray() 973 974 @*/ 975 PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr) 976 { 977 PetscErrorCode ierr; 978 PetscInt j,i,xs,ys,xm,ym,zs,zm; 979 char *iarray_start; 980 void **iptr = (void**)vptr; 981 DM_DA *dd = (DM_DA*)da->data; 982 983 PetscFunctionBegin; 984 PetscValidHeaderSpecific(da,DM_CLASSID,1); 985 if (ghosted) { 986 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 987 if (dd->arrayghostedin[i]) { 988 *iptr = dd->arrayghostedin[i]; 989 iarray_start = (char*)dd->startghostedin[i]; 990 dd->arrayghostedin[i] = PETSC_NULL; 991 dd->startghostedin[i] = PETSC_NULL; 992 993 goto done; 994 } 995 } 996 xs = dd->Xs; 997 ys = dd->Ys; 998 zs = dd->Zs; 999 xm = dd->Xe-dd->Xs; 1000 ym = dd->Ye-dd->Ys; 1001 zm = dd->Ze-dd->Zs; 1002 } else { 1003 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1004 if (dd->arrayin[i]) { 1005 *iptr = dd->arrayin[i]; 1006 iarray_start = (char*)dd->startin[i]; 1007 dd->arrayin[i] = PETSC_NULL; 1008 dd->startin[i] = PETSC_NULL; 1009 1010 goto done; 1011 } 1012 } 1013 xs = dd->xs; 1014 ys = dd->ys; 1015 zs = dd->zs; 1016 xm = dd->xe-dd->xs; 1017 ym = dd->ye-dd->ys; 1018 zm = dd->ze-dd->zs; 1019 } 1020 1021 switch (dd->dim) { 1022 case 1: { 1023 void *ptr; 1024 1025 ierr = PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1026 1027 ptr = (void*)(iarray_start - xs*sizeof(PetscScalar)); 1028 *iptr = (void*)ptr; 1029 break;} 1030 case 2: { 1031 void **ptr; 1032 1033 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1034 1035 ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*)); 1036 for(j=ys;j<ys+ym;j++) { 1037 ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs); 1038 } 1039 *iptr = (void*)ptr; 1040 break;} 1041 case 3: { 1042 void ***ptr,**bptr; 1043 1044 ierr = PetscMalloc((zm+1)*sizeof(void **)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1045 1046 ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*)); 1047 bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**)); 1048 for(i=zs;i<zs+zm;i++) { 1049 ptr[i] = bptr + ((i-zs)*ym - ys); 1050 } 1051 for (i=zs; i<zs+zm; i++) { 1052 for (j=ys; j<ys+ym; j++) { 1053 ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs); 1054 } 1055 } 1056 1057 *iptr = (void*)ptr; 1058 break;} 1059 default: 1060 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 1061 } 1062 1063 done: 1064 /* add arrays to the checked out list */ 1065 if (ghosted) { 1066 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1067 if (!dd->arrayghostedout[i]) { 1068 dd->arrayghostedout[i] = *iptr ; 1069 dd->startghostedout[i] = iarray_start; 1070 break; 1071 } 1072 } 1073 } else { 1074 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1075 if (!dd->arrayout[i]) { 1076 dd->arrayout[i] = *iptr ; 1077 dd->startout[i] = iarray_start; 1078 break; 1079 } 1080 } 1081 } 1082 PetscFunctionReturn(0); 1083 } 1084 1085 #undef __FUNCT__ 1086 #define __FUNCT__ "DMDARestoreArray" 1087 /*@C 1088 DMDARestoreArray - Restores an array of derivative types for a DMDA 1089 1090 Input Parameter: 1091 + da - information about my local patch 1092 . ghosted - do you want arrays for the ghosted or nonghosted patch 1093 - vptr - array data structured to be passed to ad_FormFunctionLocal() 1094 1095 Level: advanced 1096 1097 .seealso: DMDAGetArray(), DMDAGetAdicArray() 1098 1099 @*/ 1100 PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr) 1101 { 1102 PetscInt i; 1103 void **iptr = (void**)vptr,*iarray_start = 0; 1104 DM_DA *dd = (DM_DA*)da->data; 1105 1106 PetscFunctionBegin; 1107 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1108 if (ghosted) { 1109 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1110 if (dd->arrayghostedout[i] == *iptr) { 1111 iarray_start = dd->startghostedout[i]; 1112 dd->arrayghostedout[i] = PETSC_NULL; 1113 dd->startghostedout[i] = PETSC_NULL; 1114 break; 1115 } 1116 } 1117 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1118 if (!dd->arrayghostedin[i]){ 1119 dd->arrayghostedin[i] = *iptr; 1120 dd->startghostedin[i] = iarray_start; 1121 break; 1122 } 1123 } 1124 } else { 1125 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1126 if (dd->arrayout[i] == *iptr) { 1127 iarray_start = dd->startout[i]; 1128 dd->arrayout[i] = PETSC_NULL; 1129 dd->startout[i] = PETSC_NULL; 1130 break; 1131 } 1132 } 1133 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1134 if (!dd->arrayin[i]){ 1135 dd->arrayin[i] = *iptr; 1136 dd->startin[i] = iarray_start; 1137 break; 1138 } 1139 } 1140 } 1141 PetscFunctionReturn(0); 1142 } 1143 1144 #undef __FUNCT__ 1145 #define __FUNCT__ "DMDAGetAdicMFArray" 1146 /*@C 1147 DMDAGetAdicMFArray - Gets an array of derivative types for a DMDA for matrix-free ADIC. 1148 1149 Input Parameter: 1150 + da - information about my local patch 1151 - ghosted - do you want arrays for the ghosted or nonghosted patch? 1152 1153 Output Parameters: 1154 + vptr - array data structured to be passed to ad_FormFunctionLocal() 1155 . array_start - actual start of 1d array of all values that adiC can access directly (may be null) 1156 - tdof - total number of degrees of freedom represented in array_start (may be null) 1157 1158 Notes: 1159 The vector values are NOT initialized and may have garbage in them, so you may need 1160 to zero them. 1161 1162 This routine returns the same type of object as the DMDAVecGetArray(), except its 1163 elements are derivative types instead of PetscScalars. 1164 1165 Level: advanced 1166 1167 .seealso: DMDARestoreAdicMFArray(), DMDAGetArray(), DMDAGetAdicArray() 1168 1169 @*/ 1170 PetscErrorCode DMDAGetAdicMFArray(DM da,PetscBool ghosted,void *vptr,void *array_start,PetscInt *tdof) 1171 { 1172 PetscErrorCode ierr; 1173 PetscInt j,i,xs,ys,xm,ym,zs,zm,itdof = 0; 1174 char *iarray_start; 1175 void **iptr = (void**)vptr; 1176 DM_DA *dd = (DM_DA*)da->data; 1177 1178 PetscFunctionBegin; 1179 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1180 if (ghosted) { 1181 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1182 if (dd->admfarrayghostedin[i]) { 1183 *iptr = dd->admfarrayghostedin[i]; 1184 iarray_start = (char*)dd->admfstartghostedin[i]; 1185 itdof = dd->ghostedtdof; 1186 dd->admfarrayghostedin[i] = PETSC_NULL; 1187 dd->admfstartghostedin[i] = PETSC_NULL; 1188 1189 goto done; 1190 } 1191 } 1192 xs = dd->Xs; 1193 ys = dd->Ys; 1194 zs = dd->Zs; 1195 xm = dd->Xe-dd->Xs; 1196 ym = dd->Ye-dd->Ys; 1197 zm = dd->Ze-dd->Zs; 1198 } else { 1199 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1200 if (dd->admfarrayin[i]) { 1201 *iptr = dd->admfarrayin[i]; 1202 iarray_start = (char*)dd->admfstartin[i]; 1203 itdof = dd->tdof; 1204 dd->admfarrayin[i] = PETSC_NULL; 1205 dd->admfstartin[i] = PETSC_NULL; 1206 1207 goto done; 1208 } 1209 } 1210 xs = dd->xs; 1211 ys = dd->ys; 1212 zs = dd->zs; 1213 xm = dd->xe-dd->xs; 1214 ym = dd->ye-dd->ys; 1215 zm = dd->ze-dd->zs; 1216 } 1217 1218 switch (dd->dim) { 1219 case 1: { 1220 void *ptr; 1221 itdof = xm; 1222 1223 ierr = PetscMalloc(xm*2*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1224 1225 ptr = (void*)(iarray_start - xs*2*sizeof(PetscScalar)); 1226 *iptr = (void*)ptr; 1227 break;} 1228 case 2: { 1229 void **ptr; 1230 itdof = xm*ym; 1231 1232 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*2*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1233 1234 ptr = (void**)(iarray_start + xm*ym*2*sizeof(PetscScalar) - ys*sizeof(void*)); 1235 for(j=ys;j<ys+ym;j++) { 1236 ptr[j] = iarray_start + 2*sizeof(PetscScalar)*(xm*(j-ys) - xs); 1237 } 1238 *iptr = (void*)ptr; 1239 break;} 1240 case 3: { 1241 void ***ptr,**bptr; 1242 itdof = xm*ym*zm; 1243 1244 ierr = PetscMalloc((zm+1)*sizeof(void **)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*2*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1245 1246 ptr = (void***)(iarray_start + xm*ym*zm*2*sizeof(PetscScalar) - zs*sizeof(void*)); 1247 bptr = (void**)(iarray_start + xm*ym*zm*2*sizeof(PetscScalar) + zm*sizeof(void**)); 1248 for(i=zs;i<zs+zm;i++) { 1249 ptr[i] = bptr + ((i-zs)*ym* - ys)*sizeof(void*); 1250 } 1251 for (i=zs; i<zs+zm; i++) { 1252 for (j=ys; j<ys+ym; j++) { 1253 ptr[i][j] = iarray_start + 2*sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs); 1254 } 1255 } 1256 1257 *iptr = (void*)ptr; 1258 break;} 1259 default: 1260 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 1261 } 1262 1263 done: 1264 /* add arrays to the checked out list */ 1265 if (ghosted) { 1266 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1267 if (!dd->admfarrayghostedout[i]) { 1268 dd->admfarrayghostedout[i] = *iptr ; 1269 dd->admfstartghostedout[i] = iarray_start; 1270 dd->ghostedtdof = itdof; 1271 break; 1272 } 1273 } 1274 } else { 1275 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1276 if (!dd->admfarrayout[i]) { 1277 dd->admfarrayout[i] = *iptr ; 1278 dd->admfstartout[i] = iarray_start; 1279 dd->tdof = itdof; 1280 break; 1281 } 1282 } 1283 } 1284 if (i == DMDA_MAX_AD_ARRAYS+1) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Too many DMDA ADIC arrays obtained"); 1285 if (tdof) *tdof = itdof; 1286 if (array_start) *(void**)array_start = iarray_start; 1287 PetscFunctionReturn(0); 1288 } 1289 1290 #undef __FUNCT__ 1291 #define __FUNCT__ "DMDAGetAdicMFArray4" 1292 PetscErrorCode DMDAGetAdicMFArray4(DM da,PetscBool ghosted,void *vptr,void *array_start,PetscInt *tdof) 1293 { 1294 PetscErrorCode ierr; 1295 PetscInt j,i,xs,ys,xm,ym,itdof = 0; 1296 char *iarray_start; 1297 void **iptr = (void**)vptr; 1298 DM_DA *dd = (DM_DA*)da->data; 1299 1300 PetscFunctionBegin; 1301 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1302 if (ghosted) { 1303 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1304 if (dd->admfarrayghostedin[i]) { 1305 *iptr = dd->admfarrayghostedin[i]; 1306 iarray_start = (char*)dd->admfstartghostedin[i]; 1307 itdof = dd->ghostedtdof; 1308 dd->admfarrayghostedin[i] = PETSC_NULL; 1309 dd->admfstartghostedin[i] = PETSC_NULL; 1310 1311 goto done; 1312 } 1313 } 1314 xs = dd->Xs; 1315 ys = dd->Ys; 1316 xm = dd->Xe-dd->Xs; 1317 ym = dd->Ye-dd->Ys; 1318 } else { 1319 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1320 if (dd->admfarrayin[i]) { 1321 *iptr = dd->admfarrayin[i]; 1322 iarray_start = (char*)dd->admfstartin[i]; 1323 itdof = dd->tdof; 1324 dd->admfarrayin[i] = PETSC_NULL; 1325 dd->admfstartin[i] = PETSC_NULL; 1326 1327 goto done; 1328 } 1329 } 1330 xs = dd->xs; 1331 ys = dd->ys; 1332 xm = dd->xe-dd->xs; 1333 ym = dd->ye-dd->ys; 1334 } 1335 1336 switch (dd->dim) { 1337 case 2: { 1338 void **ptr; 1339 itdof = xm*ym; 1340 1341 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*5*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1342 1343 ptr = (void**)(iarray_start + xm*ym*5*sizeof(PetscScalar) - ys*sizeof(void*)); 1344 for(j=ys;j<ys+ym;j++) { 1345 ptr[j] = iarray_start + 5*sizeof(PetscScalar)*(xm*(j-ys) - xs); 1346 } 1347 *iptr = (void*)ptr; 1348 break;} 1349 default: 1350 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 1351 } 1352 1353 done: 1354 /* add arrays to the checked out list */ 1355 if (ghosted) { 1356 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1357 if (!dd->admfarrayghostedout[i]) { 1358 dd->admfarrayghostedout[i] = *iptr ; 1359 dd->admfstartghostedout[i] = iarray_start; 1360 dd->ghostedtdof = itdof; 1361 break; 1362 } 1363 } 1364 } else { 1365 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1366 if (!dd->admfarrayout[i]) { 1367 dd->admfarrayout[i] = *iptr ; 1368 dd->admfstartout[i] = iarray_start; 1369 dd->tdof = itdof; 1370 break; 1371 } 1372 } 1373 } 1374 if (i == DMDA_MAX_AD_ARRAYS+1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Too many DMDA ADIC arrays obtained"); 1375 if (tdof) *tdof = itdof; 1376 if (array_start) *(void**)array_start = iarray_start; 1377 PetscFunctionReturn(0); 1378 } 1379 1380 #undef __FUNCT__ 1381 #define __FUNCT__ "DMDAGetAdicMFArray9" 1382 PetscErrorCode DMDAGetAdicMFArray9(DM da,PetscBool ghosted,void *vptr,void *array_start,PetscInt *tdof) 1383 { 1384 PetscErrorCode ierr; 1385 PetscInt j,i,xs,ys,xm,ym,itdof = 0; 1386 char *iarray_start; 1387 void **iptr = (void**)vptr; 1388 DM_DA *dd = (DM_DA*)da->data; 1389 1390 PetscFunctionBegin; 1391 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1392 if (ghosted) { 1393 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1394 if (dd->admfarrayghostedin[i]) { 1395 *iptr = dd->admfarrayghostedin[i]; 1396 iarray_start = (char*)dd->admfstartghostedin[i]; 1397 itdof = dd->ghostedtdof; 1398 dd->admfarrayghostedin[i] = PETSC_NULL; 1399 dd->admfstartghostedin[i] = PETSC_NULL; 1400 1401 goto done; 1402 } 1403 } 1404 xs = dd->Xs; 1405 ys = dd->Ys; 1406 xm = dd->Xe-dd->Xs; 1407 ym = dd->Ye-dd->Ys; 1408 } else { 1409 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1410 if (dd->admfarrayin[i]) { 1411 *iptr = dd->admfarrayin[i]; 1412 iarray_start = (char*)dd->admfstartin[i]; 1413 itdof = dd->tdof; 1414 dd->admfarrayin[i] = PETSC_NULL; 1415 dd->admfstartin[i] = PETSC_NULL; 1416 1417 goto done; 1418 } 1419 } 1420 xs = dd->xs; 1421 ys = dd->ys; 1422 xm = dd->xe-dd->xs; 1423 ym = dd->ye-dd->ys; 1424 } 1425 1426 switch (dd->dim) { 1427 case 2: { 1428 void **ptr; 1429 itdof = xm*ym; 1430 1431 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*10*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1432 1433 ptr = (void**)(iarray_start + xm*ym*10*sizeof(PetscScalar) - ys*sizeof(void*)); 1434 for(j=ys;j<ys+ym;j++) { 1435 ptr[j] = iarray_start + 10*sizeof(PetscScalar)*(xm*(j-ys) - xs); 1436 } 1437 *iptr = (void*)ptr; 1438 break;} 1439 default: 1440 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 1441 } 1442 1443 done: 1444 /* add arrays to the checked out list */ 1445 if (ghosted) { 1446 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1447 if (!dd->admfarrayghostedout[i]) { 1448 dd->admfarrayghostedout[i] = *iptr ; 1449 dd->admfstartghostedout[i] = iarray_start; 1450 dd->ghostedtdof = itdof; 1451 break; 1452 } 1453 } 1454 } else { 1455 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1456 if (!dd->admfarrayout[i]) { 1457 dd->admfarrayout[i] = *iptr ; 1458 dd->admfstartout[i] = iarray_start; 1459 dd->tdof = itdof; 1460 break; 1461 } 1462 } 1463 } 1464 if (i == DMDA_MAX_AD_ARRAYS+1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Too many DMDA ADIC arrays obtained"); 1465 if (tdof) *tdof = itdof; 1466 if (array_start) *(void**)array_start = iarray_start; 1467 PetscFunctionReturn(0); 1468 } 1469 1470 #undef __FUNCT__ 1471 #define __FUNCT__ "DMDAGetAdicMFArrayb" 1472 /*@C 1473 DMDAGetAdicMFArrayb - Gets an array of derivative types for a DMDA for matrix-free ADIC. 1474 1475 Input Parameter: 1476 + da - information about my local patch 1477 - ghosted - do you want arrays for the ghosted or nonghosted patch? 1478 1479 Output Parameters: 1480 + vptr - array data structured to be passed to ad_FormFunctionLocal() 1481 . array_start - actual start of 1d array of all values that adiC can access directly (may be null) 1482 - tdof - total number of degrees of freedom represented in array_start (may be null) 1483 1484 Notes: 1485 The vector values are NOT initialized and may have garbage in them, so you may need 1486 to zero them. 1487 1488 This routine returns the same type of object as the DMDAVecGetArray(), except its 1489 elements are derivative types instead of PetscScalars. 1490 1491 Level: advanced 1492 1493 .seealso: DMDARestoreAdicMFArray(), DMDAGetArray(), DMDAGetAdicArray() 1494 1495 @*/ 1496 PetscErrorCode DMDAGetAdicMFArrayb(DM da,PetscBool ghosted,void *vptr,void *array_start,PetscInt *tdof) 1497 { 1498 PetscErrorCode ierr; 1499 PetscInt j,i,xs,ys,xm,ym,zs,zm,itdof = 0; 1500 char *iarray_start; 1501 void **iptr = (void**)vptr; 1502 DM_DA *dd = (DM_DA*)da->data; 1503 PetscInt bs = dd->w,bs1 = bs+1; 1504 1505 PetscFunctionBegin; 1506 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1507 if (ghosted) { 1508 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1509 if (dd->admfarrayghostedin[i]) { 1510 *iptr = dd->admfarrayghostedin[i]; 1511 iarray_start = (char*)dd->admfstartghostedin[i]; 1512 itdof = dd->ghostedtdof; 1513 dd->admfarrayghostedin[i] = PETSC_NULL; 1514 dd->admfstartghostedin[i] = PETSC_NULL; 1515 1516 goto done; 1517 } 1518 } 1519 xs = dd->Xs; 1520 ys = dd->Ys; 1521 zs = dd->Zs; 1522 xm = dd->Xe-dd->Xs; 1523 ym = dd->Ye-dd->Ys; 1524 zm = dd->Ze-dd->Zs; 1525 } else { 1526 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1527 if (dd->admfarrayin[i]) { 1528 *iptr = dd->admfarrayin[i]; 1529 iarray_start = (char*)dd->admfstartin[i]; 1530 itdof = dd->tdof; 1531 dd->admfarrayin[i] = PETSC_NULL; 1532 dd->admfstartin[i] = PETSC_NULL; 1533 1534 goto done; 1535 } 1536 } 1537 xs = dd->xs; 1538 ys = dd->ys; 1539 zs = dd->zs; 1540 xm = dd->xe-dd->xs; 1541 ym = dd->ye-dd->ys; 1542 zm = dd->ze-dd->zs; 1543 } 1544 1545 switch (dd->dim) { 1546 case 1: { 1547 void *ptr; 1548 itdof = xm; 1549 1550 ierr = PetscMalloc(xm*bs1*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1551 1552 ptr = (void*)(iarray_start - xs*bs1*sizeof(PetscScalar)); 1553 *iptr = (void*)ptr; 1554 break;} 1555 case 2: { 1556 void **ptr; 1557 itdof = xm*ym; 1558 1559 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*bs1*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1560 1561 ptr = (void**)(iarray_start + xm*ym*bs1*sizeof(PetscScalar) - ys*sizeof(void*)); 1562 for(j=ys;j<ys+ym;j++) { 1563 ptr[j] = iarray_start + bs1*sizeof(PetscScalar)*(xm*(j-ys) - xs); 1564 } 1565 *iptr = (void*)ptr; 1566 break;} 1567 case 3: { 1568 void ***ptr,**bptr; 1569 itdof = xm*ym*zm; 1570 1571 ierr = PetscMalloc((zm+1)*sizeof(void **)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*bs1*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1572 1573 ptr = (void***)(iarray_start + xm*ym*zm*2*sizeof(PetscScalar) - zs*sizeof(void*)); 1574 bptr = (void**)(iarray_start + xm*ym*zm*2*sizeof(PetscScalar) + zm*sizeof(void**)); 1575 for(i=zs;i<zs+zm;i++) { 1576 ptr[i] = bptr + ((i-zs)*ym* - ys)*sizeof(void*); 1577 } 1578 for (i=zs; i<zs+zm; i++) { 1579 for (j=ys; j<ys+ym; j++) { 1580 ptr[i][j] = iarray_start + bs1*sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs); 1581 } 1582 } 1583 1584 *iptr = (void*)ptr; 1585 break;} 1586 default: 1587 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 1588 } 1589 1590 done: 1591 /* add arrays to the checked out list */ 1592 if (ghosted) { 1593 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1594 if (!dd->admfarrayghostedout[i]) { 1595 dd->admfarrayghostedout[i] = *iptr ; 1596 dd->admfstartghostedout[i] = iarray_start; 1597 dd->ghostedtdof = itdof; 1598 break; 1599 } 1600 } 1601 } else { 1602 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1603 if (!dd->admfarrayout[i]) { 1604 dd->admfarrayout[i] = *iptr ; 1605 dd->admfstartout[i] = iarray_start; 1606 dd->tdof = itdof; 1607 break; 1608 } 1609 } 1610 } 1611 if (i == DMDA_MAX_AD_ARRAYS+1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Too many DMDA ADIC arrays obtained"); 1612 if (tdof) *tdof = itdof; 1613 if (array_start) *(void**)array_start = iarray_start; 1614 PetscFunctionReturn(0); 1615 } 1616 1617 #undef __FUNCT__ 1618 #define __FUNCT__ "DMDARestoreAdicMFArray" 1619 /*@C 1620 DMDARestoreAdicMFArray - Restores an array of derivative types for a DMDA. 1621 1622 Input Parameter: 1623 + da - information about my local patch 1624 - ghosted - do you want arrays for the ghosted or nonghosted patch? 1625 1626 Output Parameters: 1627 + ptr - array data structure to be passed to ad_FormFunctionLocal() 1628 . array_start - actual start of 1d array of all values that adiC can access directly 1629 - tdof - total number of degrees of freedom represented in array_start 1630 1631 Level: advanced 1632 1633 .seealso: DMDAGetAdicArray() 1634 1635 @*/ 1636 PetscErrorCode DMDARestoreAdicMFArray(DM da,PetscBool ghosted,void *vptr,void *array_start,PetscInt *tdof) 1637 { 1638 PetscInt i; 1639 void **iptr = (void**)vptr,*iarray_start = 0; 1640 DM_DA *dd = (DM_DA*)da->data; 1641 1642 PetscFunctionBegin; 1643 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1644 if (ghosted) { 1645 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1646 if (dd->admfarrayghostedout[i] == *iptr) { 1647 iarray_start = dd->admfstartghostedout[i]; 1648 dd->admfarrayghostedout[i] = PETSC_NULL; 1649 dd->admfstartghostedout[i] = PETSC_NULL; 1650 break; 1651 } 1652 } 1653 if (!iarray_start) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Could not find array in checkout list"); 1654 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1655 if (!dd->admfarrayghostedin[i]){ 1656 dd->admfarrayghostedin[i] = *iptr; 1657 dd->admfstartghostedin[i] = iarray_start; 1658 break; 1659 } 1660 } 1661 } else { 1662 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1663 if (dd->admfarrayout[i] == *iptr) { 1664 iarray_start = dd->admfstartout[i]; 1665 dd->admfarrayout[i] = PETSC_NULL; 1666 dd->admfstartout[i] = PETSC_NULL; 1667 break; 1668 } 1669 } 1670 if (!iarray_start) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Could not find array in checkout list"); 1671 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1672 if (!dd->admfarrayin[i]){ 1673 dd->admfarrayin[i] = *iptr; 1674 dd->admfstartin[i] = iarray_start; 1675 break; 1676 } 1677 } 1678 } 1679 PetscFunctionReturn(0); 1680 } 1681 1682 #undef __FUNCT__ 1683 #define __FUNCT__ "admf_DAGetArray" 1684 PetscErrorCode admf_DAGetArray(DM da,PetscBool ghosted,void *iptr) 1685 { 1686 PetscErrorCode ierr; 1687 PetscFunctionBegin; 1688 ierr = DMDAGetAdicMFArray(da,ghosted,iptr,0,0);CHKERRQ(ierr); 1689 PetscFunctionReturn(0); 1690 } 1691 1692 #undef __FUNCT__ 1693 #define __FUNCT__ "admf_DARestoreArray" 1694 PetscErrorCode admf_DARestoreArray(DM da,PetscBool ghosted,void *iptr) 1695 { 1696 PetscErrorCode ierr; 1697 PetscFunctionBegin; 1698 ierr = DMDARestoreAdicMFArray(da,ghosted,iptr,0,0);CHKERRQ(ierr); 1699 PetscFunctionReturn(0); 1700 } 1701 1702