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 /* TODO: Remove duplication in leaf determination */ 704 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); 705 ierr = PetscSFCreate(((PetscObject) dm)->comm, &sf);CHKERRQ(ierr); 706 ierr = PetscSFSetGraph(sf, pEnd, nleaves, localPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 707 /* Create global section */ 708 ierr = PetscSectionCreateGlobalSection(dm->defaultSection, sf, &dm->defaultGlobalSection);CHKERRQ(ierr); 709 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 710 /* Create default SF */ 711 ierr = DMCreateDefaultSF(dm, dm->defaultSection, dm->defaultGlobalSection);CHKERRQ(ierr); 712 PetscFunctionReturn(0); 713 } 714 715 /* ------------------------------------------------------------------- */ 716 #if defined(PETSC_HAVE_ADIC) 717 718 EXTERN_C_BEGIN 719 #include <adic/ad_utils.h> 720 EXTERN_C_END 721 722 #undef __FUNCT__ 723 #define __FUNCT__ "DMDAGetAdicArray" 724 /*@C 725 DMDAGetAdicArray - Gets an array of derivative types for a DMDA 726 727 Input Parameter: 728 + da - information about my local patch 729 - ghosted - do you want arrays for the ghosted or nonghosted patch 730 731 Output Parameters: 732 + vptr - array data structured to be passed to ad_FormFunctionLocal() 733 . array_start - actual start of 1d array of all values that adiC can access directly (may be null) 734 - tdof - total number of degrees of freedom represented in array_start (may be null) 735 736 Notes: 737 The vector values are NOT initialized and may have garbage in them, so you may need 738 to zero them. 739 740 Returns the same type of object as the DMDAVecGetArray() except its elements are 741 derivative types instead of PetscScalars 742 743 Level: advanced 744 745 .seealso: DMDARestoreAdicArray() 746 747 @*/ 748 PetscErrorCode DMDAGetAdicArray(DM da,PetscBool ghosted,void *vptr,void *array_start,PetscInt *tdof) 749 { 750 PetscErrorCode ierr; 751 PetscInt j,i,deriv_type_size,xs,ys,xm,ym,zs,zm,itdof; 752 char *iarray_start; 753 void **iptr = (void**)vptr; 754 DM_DA *dd = (DM_DA*)da->data; 755 756 PetscFunctionBegin; 757 PetscValidHeaderSpecific(da,DM_CLASSID,1); 758 if (ghosted) { 759 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 760 if (dd->adarrayghostedin[i]) { 761 *iptr = dd->adarrayghostedin[i]; 762 iarray_start = (char*)dd->adstartghostedin[i]; 763 itdof = dd->ghostedtdof; 764 dd->adarrayghostedin[i] = PETSC_NULL; 765 dd->adstartghostedin[i] = PETSC_NULL; 766 767 goto done; 768 } 769 } 770 xs = dd->Xs; 771 ys = dd->Ys; 772 zs = dd->Zs; 773 xm = dd->Xe-dd->Xs; 774 ym = dd->Ye-dd->Ys; 775 zm = dd->Ze-dd->Zs; 776 } else { 777 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 778 if (dd->adarrayin[i]) { 779 *iptr = dd->adarrayin[i]; 780 iarray_start = (char*)dd->adstartin[i]; 781 itdof = dd->tdof; 782 dd->adarrayin[i] = PETSC_NULL; 783 dd->adstartin[i] = PETSC_NULL; 784 785 goto done; 786 } 787 } 788 xs = dd->xs; 789 ys = dd->ys; 790 zs = dd->zs; 791 xm = dd->xe-dd->xs; 792 ym = dd->ye-dd->ys; 793 zm = dd->ze-dd->zs; 794 } 795 deriv_type_size = PetscADGetDerivTypeSize(); 796 797 switch (dd->dim) { 798 case 1: { 799 void *ptr; 800 itdof = xm; 801 802 ierr = PetscMalloc(xm*deriv_type_size,&iarray_start);CHKERRQ(ierr); 803 804 ptr = (void*)(iarray_start - xs*deriv_type_size); 805 *iptr = (void*)ptr; 806 break;} 807 case 2: { 808 void **ptr; 809 itdof = xm*ym; 810 811 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*deriv_type_size,&iarray_start);CHKERRQ(ierr); 812 813 ptr = (void**)(iarray_start + xm*ym*deriv_type_size - ys*sizeof(void*)); 814 for(j=ys;j<ys+ym;j++) { 815 ptr[j] = iarray_start + deriv_type_size*(xm*(j-ys) - xs); 816 } 817 *iptr = (void*)ptr; 818 break;} 819 case 3: { 820 void ***ptr,**bptr; 821 itdof = xm*ym*zm; 822 823 ierr = PetscMalloc((zm+1)*sizeof(void **)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*deriv_type_size,&iarray_start);CHKERRQ(ierr); 824 825 ptr = (void***)(iarray_start + xm*ym*zm*deriv_type_size - zs*sizeof(void*)); 826 bptr = (void**)(iarray_start + xm*ym*zm*deriv_type_size + zm*sizeof(void**)); 827 828 for(i=zs;i<zs+zm;i++) { 829 ptr[i] = bptr + ((i-zs)*ym - ys); 830 } 831 for (i=zs; i<zs+zm; i++) { 832 for (j=ys; j<ys+ym; j++) { 833 ptr[i][j] = iarray_start + deriv_type_size*(xm*ym*(i-zs) + xm*(j-ys) - xs); 834 } 835 } 836 837 *iptr = (void*)ptr; 838 break;} 839 default: 840 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 841 } 842 843 done: 844 /* add arrays to the checked out list */ 845 if (ghosted) { 846 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 847 if (!dd->adarrayghostedout[i]) { 848 dd->adarrayghostedout[i] = *iptr ; 849 dd->adstartghostedout[i] = iarray_start; 850 dd->ghostedtdof = itdof; 851 break; 852 } 853 } 854 } else { 855 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 856 if (!dd->adarrayout[i]) { 857 dd->adarrayout[i] = *iptr ; 858 dd->adstartout[i] = iarray_start; 859 dd->tdof = itdof; 860 break; 861 } 862 } 863 } 864 if (i == DMDA_MAX_AD_ARRAYS+1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Too many DMDA ADIC arrays obtained"); 865 if (tdof) *tdof = itdof; 866 if (array_start) *(void**)array_start = iarray_start; 867 PetscFunctionReturn(0); 868 } 869 870 #undef __FUNCT__ 871 #define __FUNCT__ "DMDARestoreAdicArray" 872 /*@C 873 DMDARestoreAdicArray - Restores an array of derivative types for a DMDA 874 875 Input Parameter: 876 + da - information about my local patch 877 - ghosted - do you want arrays for the ghosted or nonghosted patch 878 879 Output Parameters: 880 + ptr - array data structured to be passed to ad_FormFunctionLocal() 881 . array_start - actual start of 1d array of all values that adiC can access directly 882 - tdof - total number of degrees of freedom represented in array_start 883 884 Level: advanced 885 886 .seealso: DMDAGetAdicArray() 887 888 @*/ 889 PetscErrorCode DMDARestoreAdicArray(DM da,PetscBool ghosted,void *ptr,void *array_start,PetscInt *tdof) 890 { 891 PetscInt i; 892 void **iptr = (void**)ptr,iarray_start = 0; 893 894 PetscFunctionBegin; 895 PetscValidHeaderSpecific(da,DM_CLASSID,1); 896 if (ghosted) { 897 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 898 if (dd->adarrayghostedout[i] == *iptr) { 899 iarray_start = dd->adstartghostedout[i]; 900 dd->adarrayghostedout[i] = PETSC_NULL; 901 dd->adstartghostedout[i] = PETSC_NULL; 902 break; 903 } 904 } 905 if (!iarray_start) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Could not find array in checkout list"); 906 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 907 if (!dd->adarrayghostedin[i]){ 908 dd->adarrayghostedin[i] = *iptr; 909 dd->adstartghostedin[i] = iarray_start; 910 break; 911 } 912 } 913 } else { 914 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 915 if (dd->adarrayout[i] == *iptr) { 916 iarray_start = dd->adstartout[i]; 917 dd->adarrayout[i] = PETSC_NULL; 918 dd->adstartout[i] = PETSC_NULL; 919 break; 920 } 921 } 922 if (!iarray_start) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Could not find array in checkout list"); 923 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 924 if (!dd->adarrayin[i]){ 925 dd->adarrayin[i] = *iptr; 926 dd->adstartin[i] = iarray_start; 927 break; 928 } 929 } 930 } 931 PetscFunctionReturn(0); 932 } 933 934 #undef __FUNCT__ 935 #define __FUNCT__ "ad_DAGetArray" 936 PetscErrorCode ad_DAGetArray(DM da,PetscBool ghosted,void *iptr) 937 { 938 PetscErrorCode ierr; 939 PetscFunctionBegin; 940 ierr = DMDAGetAdicArray(da,ghosted,iptr,0,0);CHKERRQ(ierr); 941 PetscFunctionReturn(0); 942 } 943 944 #undef __FUNCT__ 945 #define __FUNCT__ "ad_DARestoreArray" 946 PetscErrorCode ad_DARestoreArray(DM da,PetscBool ghosted,void *iptr) 947 { 948 PetscErrorCode ierr; 949 PetscFunctionBegin; 950 ierr = DMDARestoreAdicArray(da,ghosted,iptr,0,0);CHKERRQ(ierr); 951 PetscFunctionReturn(0); 952 } 953 954 #endif 955 956 #undef __FUNCT__ 957 #define __FUNCT__ "DMDAGetArray" 958 /*@C 959 DMDAGetArray - Gets a work array for a DMDA 960 961 Input Parameter: 962 + da - information about my local patch 963 - ghosted - do you want arrays for the ghosted or nonghosted patch 964 965 Output Parameters: 966 . vptr - array data structured 967 968 Note: The vector values are NOT initialized and may have garbage in them, so you may need 969 to zero them. 970 971 Level: advanced 972 973 .seealso: DMDARestoreArray(), DMDAGetAdicArray() 974 975 @*/ 976 PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr) 977 { 978 PetscErrorCode ierr; 979 PetscInt j,i,xs,ys,xm,ym,zs,zm; 980 char *iarray_start; 981 void **iptr = (void**)vptr; 982 DM_DA *dd = (DM_DA*)da->data; 983 984 PetscFunctionBegin; 985 PetscValidHeaderSpecific(da,DM_CLASSID,1); 986 if (ghosted) { 987 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 988 if (dd->arrayghostedin[i]) { 989 *iptr = dd->arrayghostedin[i]; 990 iarray_start = (char*)dd->startghostedin[i]; 991 dd->arrayghostedin[i] = PETSC_NULL; 992 dd->startghostedin[i] = PETSC_NULL; 993 994 goto done; 995 } 996 } 997 xs = dd->Xs; 998 ys = dd->Ys; 999 zs = dd->Zs; 1000 xm = dd->Xe-dd->Xs; 1001 ym = dd->Ye-dd->Ys; 1002 zm = dd->Ze-dd->Zs; 1003 } else { 1004 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1005 if (dd->arrayin[i]) { 1006 *iptr = dd->arrayin[i]; 1007 iarray_start = (char*)dd->startin[i]; 1008 dd->arrayin[i] = PETSC_NULL; 1009 dd->startin[i] = PETSC_NULL; 1010 1011 goto done; 1012 } 1013 } 1014 xs = dd->xs; 1015 ys = dd->ys; 1016 zs = dd->zs; 1017 xm = dd->xe-dd->xs; 1018 ym = dd->ye-dd->ys; 1019 zm = dd->ze-dd->zs; 1020 } 1021 1022 switch (dd->dim) { 1023 case 1: { 1024 void *ptr; 1025 1026 ierr = PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1027 1028 ptr = (void*)(iarray_start - xs*sizeof(PetscScalar)); 1029 *iptr = (void*)ptr; 1030 break;} 1031 case 2: { 1032 void **ptr; 1033 1034 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1035 1036 ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*)); 1037 for(j=ys;j<ys+ym;j++) { 1038 ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs); 1039 } 1040 *iptr = (void*)ptr; 1041 break;} 1042 case 3: { 1043 void ***ptr,**bptr; 1044 1045 ierr = PetscMalloc((zm+1)*sizeof(void **)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1046 1047 ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*)); 1048 bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**)); 1049 for(i=zs;i<zs+zm;i++) { 1050 ptr[i] = bptr + ((i-zs)*ym - ys); 1051 } 1052 for (i=zs; i<zs+zm; i++) { 1053 for (j=ys; j<ys+ym; j++) { 1054 ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs); 1055 } 1056 } 1057 1058 *iptr = (void*)ptr; 1059 break;} 1060 default: 1061 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 1062 } 1063 1064 done: 1065 /* add arrays to the checked out list */ 1066 if (ghosted) { 1067 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1068 if (!dd->arrayghostedout[i]) { 1069 dd->arrayghostedout[i] = *iptr ; 1070 dd->startghostedout[i] = iarray_start; 1071 break; 1072 } 1073 } 1074 } else { 1075 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1076 if (!dd->arrayout[i]) { 1077 dd->arrayout[i] = *iptr ; 1078 dd->startout[i] = iarray_start; 1079 break; 1080 } 1081 } 1082 } 1083 PetscFunctionReturn(0); 1084 } 1085 1086 #undef __FUNCT__ 1087 #define __FUNCT__ "DMDARestoreArray" 1088 /*@C 1089 DMDARestoreArray - Restores an array of derivative types for a DMDA 1090 1091 Input Parameter: 1092 + da - information about my local patch 1093 . ghosted - do you want arrays for the ghosted or nonghosted patch 1094 - vptr - array data structured to be passed to ad_FormFunctionLocal() 1095 1096 Level: advanced 1097 1098 .seealso: DMDAGetArray(), DMDAGetAdicArray() 1099 1100 @*/ 1101 PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr) 1102 { 1103 PetscInt i; 1104 void **iptr = (void**)vptr,*iarray_start = 0; 1105 DM_DA *dd = (DM_DA*)da->data; 1106 1107 PetscFunctionBegin; 1108 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1109 if (ghosted) { 1110 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1111 if (dd->arrayghostedout[i] == *iptr) { 1112 iarray_start = dd->startghostedout[i]; 1113 dd->arrayghostedout[i] = PETSC_NULL; 1114 dd->startghostedout[i] = PETSC_NULL; 1115 break; 1116 } 1117 } 1118 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1119 if (!dd->arrayghostedin[i]){ 1120 dd->arrayghostedin[i] = *iptr; 1121 dd->startghostedin[i] = iarray_start; 1122 break; 1123 } 1124 } 1125 } else { 1126 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1127 if (dd->arrayout[i] == *iptr) { 1128 iarray_start = dd->startout[i]; 1129 dd->arrayout[i] = PETSC_NULL; 1130 dd->startout[i] = PETSC_NULL; 1131 break; 1132 } 1133 } 1134 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1135 if (!dd->arrayin[i]){ 1136 dd->arrayin[i] = *iptr; 1137 dd->startin[i] = iarray_start; 1138 break; 1139 } 1140 } 1141 } 1142 PetscFunctionReturn(0); 1143 } 1144 1145 #undef __FUNCT__ 1146 #define __FUNCT__ "DMDAGetAdicMFArray" 1147 /*@C 1148 DMDAGetAdicMFArray - Gets an array of derivative types for a DMDA for matrix-free ADIC. 1149 1150 Input Parameter: 1151 + da - information about my local patch 1152 - ghosted - do you want arrays for the ghosted or nonghosted patch? 1153 1154 Output Parameters: 1155 + vptr - array data structured to be passed to ad_FormFunctionLocal() 1156 . array_start - actual start of 1d array of all values that adiC can access directly (may be null) 1157 - tdof - total number of degrees of freedom represented in array_start (may be null) 1158 1159 Notes: 1160 The vector values are NOT initialized and may have garbage in them, so you may need 1161 to zero them. 1162 1163 This routine returns the same type of object as the DMDAVecGetArray(), except its 1164 elements are derivative types instead of PetscScalars. 1165 1166 Level: advanced 1167 1168 .seealso: DMDARestoreAdicMFArray(), DMDAGetArray(), DMDAGetAdicArray() 1169 1170 @*/ 1171 PetscErrorCode DMDAGetAdicMFArray(DM da,PetscBool ghosted,void *vptr,void *array_start,PetscInt *tdof) 1172 { 1173 PetscErrorCode ierr; 1174 PetscInt j,i,xs,ys,xm,ym,zs,zm,itdof = 0; 1175 char *iarray_start; 1176 void **iptr = (void**)vptr; 1177 DM_DA *dd = (DM_DA*)da->data; 1178 1179 PetscFunctionBegin; 1180 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1181 if (ghosted) { 1182 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1183 if (dd->admfarrayghostedin[i]) { 1184 *iptr = dd->admfarrayghostedin[i]; 1185 iarray_start = (char*)dd->admfstartghostedin[i]; 1186 itdof = dd->ghostedtdof; 1187 dd->admfarrayghostedin[i] = PETSC_NULL; 1188 dd->admfstartghostedin[i] = PETSC_NULL; 1189 1190 goto done; 1191 } 1192 } 1193 xs = dd->Xs; 1194 ys = dd->Ys; 1195 zs = dd->Zs; 1196 xm = dd->Xe-dd->Xs; 1197 ym = dd->Ye-dd->Ys; 1198 zm = dd->Ze-dd->Zs; 1199 } else { 1200 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1201 if (dd->admfarrayin[i]) { 1202 *iptr = dd->admfarrayin[i]; 1203 iarray_start = (char*)dd->admfstartin[i]; 1204 itdof = dd->tdof; 1205 dd->admfarrayin[i] = PETSC_NULL; 1206 dd->admfstartin[i] = PETSC_NULL; 1207 1208 goto done; 1209 } 1210 } 1211 xs = dd->xs; 1212 ys = dd->ys; 1213 zs = dd->zs; 1214 xm = dd->xe-dd->xs; 1215 ym = dd->ye-dd->ys; 1216 zm = dd->ze-dd->zs; 1217 } 1218 1219 switch (dd->dim) { 1220 case 1: { 1221 void *ptr; 1222 itdof = xm; 1223 1224 ierr = PetscMalloc(xm*2*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1225 1226 ptr = (void*)(iarray_start - xs*2*sizeof(PetscScalar)); 1227 *iptr = (void*)ptr; 1228 break;} 1229 case 2: { 1230 void **ptr; 1231 itdof = xm*ym; 1232 1233 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*2*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1234 1235 ptr = (void**)(iarray_start + xm*ym*2*sizeof(PetscScalar) - ys*sizeof(void*)); 1236 for(j=ys;j<ys+ym;j++) { 1237 ptr[j] = iarray_start + 2*sizeof(PetscScalar)*(xm*(j-ys) - xs); 1238 } 1239 *iptr = (void*)ptr; 1240 break;} 1241 case 3: { 1242 void ***ptr,**bptr; 1243 itdof = xm*ym*zm; 1244 1245 ierr = PetscMalloc((zm+1)*sizeof(void **)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*2*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1246 1247 ptr = (void***)(iarray_start + xm*ym*zm*2*sizeof(PetscScalar) - zs*sizeof(void*)); 1248 bptr = (void**)(iarray_start + xm*ym*zm*2*sizeof(PetscScalar) + zm*sizeof(void**)); 1249 for(i=zs;i<zs+zm;i++) { 1250 ptr[i] = bptr + ((i-zs)*ym* - ys)*sizeof(void*); 1251 } 1252 for (i=zs; i<zs+zm; i++) { 1253 for (j=ys; j<ys+ym; j++) { 1254 ptr[i][j] = iarray_start + 2*sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs); 1255 } 1256 } 1257 1258 *iptr = (void*)ptr; 1259 break;} 1260 default: 1261 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 1262 } 1263 1264 done: 1265 /* add arrays to the checked out list */ 1266 if (ghosted) { 1267 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1268 if (!dd->admfarrayghostedout[i]) { 1269 dd->admfarrayghostedout[i] = *iptr ; 1270 dd->admfstartghostedout[i] = iarray_start; 1271 dd->ghostedtdof = itdof; 1272 break; 1273 } 1274 } 1275 } else { 1276 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1277 if (!dd->admfarrayout[i]) { 1278 dd->admfarrayout[i] = *iptr ; 1279 dd->admfstartout[i] = iarray_start; 1280 dd->tdof = itdof; 1281 break; 1282 } 1283 } 1284 } 1285 if (i == DMDA_MAX_AD_ARRAYS+1) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Too many DMDA ADIC arrays obtained"); 1286 if (tdof) *tdof = itdof; 1287 if (array_start) *(void**)array_start = iarray_start; 1288 PetscFunctionReturn(0); 1289 } 1290 1291 #undef __FUNCT__ 1292 #define __FUNCT__ "DMDAGetAdicMFArray4" 1293 PetscErrorCode DMDAGetAdicMFArray4(DM da,PetscBool ghosted,void *vptr,void *array_start,PetscInt *tdof) 1294 { 1295 PetscErrorCode ierr; 1296 PetscInt j,i,xs,ys,xm,ym,itdof = 0; 1297 char *iarray_start; 1298 void **iptr = (void**)vptr; 1299 DM_DA *dd = (DM_DA*)da->data; 1300 1301 PetscFunctionBegin; 1302 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1303 if (ghosted) { 1304 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1305 if (dd->admfarrayghostedin[i]) { 1306 *iptr = dd->admfarrayghostedin[i]; 1307 iarray_start = (char*)dd->admfstartghostedin[i]; 1308 itdof = dd->ghostedtdof; 1309 dd->admfarrayghostedin[i] = PETSC_NULL; 1310 dd->admfstartghostedin[i] = PETSC_NULL; 1311 1312 goto done; 1313 } 1314 } 1315 xs = dd->Xs; 1316 ys = dd->Ys; 1317 xm = dd->Xe-dd->Xs; 1318 ym = dd->Ye-dd->Ys; 1319 } else { 1320 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1321 if (dd->admfarrayin[i]) { 1322 *iptr = dd->admfarrayin[i]; 1323 iarray_start = (char*)dd->admfstartin[i]; 1324 itdof = dd->tdof; 1325 dd->admfarrayin[i] = PETSC_NULL; 1326 dd->admfstartin[i] = PETSC_NULL; 1327 1328 goto done; 1329 } 1330 } 1331 xs = dd->xs; 1332 ys = dd->ys; 1333 xm = dd->xe-dd->xs; 1334 ym = dd->ye-dd->ys; 1335 } 1336 1337 switch (dd->dim) { 1338 case 2: { 1339 void **ptr; 1340 itdof = xm*ym; 1341 1342 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*5*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1343 1344 ptr = (void**)(iarray_start + xm*ym*5*sizeof(PetscScalar) - ys*sizeof(void*)); 1345 for(j=ys;j<ys+ym;j++) { 1346 ptr[j] = iarray_start + 5*sizeof(PetscScalar)*(xm*(j-ys) - xs); 1347 } 1348 *iptr = (void*)ptr; 1349 break;} 1350 default: 1351 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 1352 } 1353 1354 done: 1355 /* add arrays to the checked out list */ 1356 if (ghosted) { 1357 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1358 if (!dd->admfarrayghostedout[i]) { 1359 dd->admfarrayghostedout[i] = *iptr ; 1360 dd->admfstartghostedout[i] = iarray_start; 1361 dd->ghostedtdof = itdof; 1362 break; 1363 } 1364 } 1365 } else { 1366 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1367 if (!dd->admfarrayout[i]) { 1368 dd->admfarrayout[i] = *iptr ; 1369 dd->admfstartout[i] = iarray_start; 1370 dd->tdof = itdof; 1371 break; 1372 } 1373 } 1374 } 1375 if (i == DMDA_MAX_AD_ARRAYS+1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Too many DMDA ADIC arrays obtained"); 1376 if (tdof) *tdof = itdof; 1377 if (array_start) *(void**)array_start = iarray_start; 1378 PetscFunctionReturn(0); 1379 } 1380 1381 #undef __FUNCT__ 1382 #define __FUNCT__ "DMDAGetAdicMFArray9" 1383 PetscErrorCode DMDAGetAdicMFArray9(DM da,PetscBool ghosted,void *vptr,void *array_start,PetscInt *tdof) 1384 { 1385 PetscErrorCode ierr; 1386 PetscInt j,i,xs,ys,xm,ym,itdof = 0; 1387 char *iarray_start; 1388 void **iptr = (void**)vptr; 1389 DM_DA *dd = (DM_DA*)da->data; 1390 1391 PetscFunctionBegin; 1392 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1393 if (ghosted) { 1394 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1395 if (dd->admfarrayghostedin[i]) { 1396 *iptr = dd->admfarrayghostedin[i]; 1397 iarray_start = (char*)dd->admfstartghostedin[i]; 1398 itdof = dd->ghostedtdof; 1399 dd->admfarrayghostedin[i] = PETSC_NULL; 1400 dd->admfstartghostedin[i] = PETSC_NULL; 1401 1402 goto done; 1403 } 1404 } 1405 xs = dd->Xs; 1406 ys = dd->Ys; 1407 xm = dd->Xe-dd->Xs; 1408 ym = dd->Ye-dd->Ys; 1409 } else { 1410 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1411 if (dd->admfarrayin[i]) { 1412 *iptr = dd->admfarrayin[i]; 1413 iarray_start = (char*)dd->admfstartin[i]; 1414 itdof = dd->tdof; 1415 dd->admfarrayin[i] = PETSC_NULL; 1416 dd->admfstartin[i] = PETSC_NULL; 1417 1418 goto done; 1419 } 1420 } 1421 xs = dd->xs; 1422 ys = dd->ys; 1423 xm = dd->xe-dd->xs; 1424 ym = dd->ye-dd->ys; 1425 } 1426 1427 switch (dd->dim) { 1428 case 2: { 1429 void **ptr; 1430 itdof = xm*ym; 1431 1432 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*10*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1433 1434 ptr = (void**)(iarray_start + xm*ym*10*sizeof(PetscScalar) - ys*sizeof(void*)); 1435 for(j=ys;j<ys+ym;j++) { 1436 ptr[j] = iarray_start + 10*sizeof(PetscScalar)*(xm*(j-ys) - xs); 1437 } 1438 *iptr = (void*)ptr; 1439 break;} 1440 default: 1441 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 1442 } 1443 1444 done: 1445 /* add arrays to the checked out list */ 1446 if (ghosted) { 1447 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1448 if (!dd->admfarrayghostedout[i]) { 1449 dd->admfarrayghostedout[i] = *iptr ; 1450 dd->admfstartghostedout[i] = iarray_start; 1451 dd->ghostedtdof = itdof; 1452 break; 1453 } 1454 } 1455 } else { 1456 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1457 if (!dd->admfarrayout[i]) { 1458 dd->admfarrayout[i] = *iptr ; 1459 dd->admfstartout[i] = iarray_start; 1460 dd->tdof = itdof; 1461 break; 1462 } 1463 } 1464 } 1465 if (i == DMDA_MAX_AD_ARRAYS+1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Too many DMDA ADIC arrays obtained"); 1466 if (tdof) *tdof = itdof; 1467 if (array_start) *(void**)array_start = iarray_start; 1468 PetscFunctionReturn(0); 1469 } 1470 1471 #undef __FUNCT__ 1472 #define __FUNCT__ "DMDAGetAdicMFArrayb" 1473 /*@C 1474 DMDAGetAdicMFArrayb - Gets an array of derivative types for a DMDA for matrix-free ADIC. 1475 1476 Input Parameter: 1477 + da - information about my local patch 1478 - ghosted - do you want arrays for the ghosted or nonghosted patch? 1479 1480 Output Parameters: 1481 + vptr - array data structured to be passed to ad_FormFunctionLocal() 1482 . array_start - actual start of 1d array of all values that adiC can access directly (may be null) 1483 - tdof - total number of degrees of freedom represented in array_start (may be null) 1484 1485 Notes: 1486 The vector values are NOT initialized and may have garbage in them, so you may need 1487 to zero them. 1488 1489 This routine returns the same type of object as the DMDAVecGetArray(), except its 1490 elements are derivative types instead of PetscScalars. 1491 1492 Level: advanced 1493 1494 .seealso: DMDARestoreAdicMFArray(), DMDAGetArray(), DMDAGetAdicArray() 1495 1496 @*/ 1497 PetscErrorCode DMDAGetAdicMFArrayb(DM da,PetscBool ghosted,void *vptr,void *array_start,PetscInt *tdof) 1498 { 1499 PetscErrorCode ierr; 1500 PetscInt j,i,xs,ys,xm,ym,zs,zm,itdof = 0; 1501 char *iarray_start; 1502 void **iptr = (void**)vptr; 1503 DM_DA *dd = (DM_DA*)da->data; 1504 PetscInt bs = dd->w,bs1 = bs+1; 1505 1506 PetscFunctionBegin; 1507 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1508 if (ghosted) { 1509 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1510 if (dd->admfarrayghostedin[i]) { 1511 *iptr = dd->admfarrayghostedin[i]; 1512 iarray_start = (char*)dd->admfstartghostedin[i]; 1513 itdof = dd->ghostedtdof; 1514 dd->admfarrayghostedin[i] = PETSC_NULL; 1515 dd->admfstartghostedin[i] = PETSC_NULL; 1516 1517 goto done; 1518 } 1519 } 1520 xs = dd->Xs; 1521 ys = dd->Ys; 1522 zs = dd->Zs; 1523 xm = dd->Xe-dd->Xs; 1524 ym = dd->Ye-dd->Ys; 1525 zm = dd->Ze-dd->Zs; 1526 } else { 1527 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1528 if (dd->admfarrayin[i]) { 1529 *iptr = dd->admfarrayin[i]; 1530 iarray_start = (char*)dd->admfstartin[i]; 1531 itdof = dd->tdof; 1532 dd->admfarrayin[i] = PETSC_NULL; 1533 dd->admfstartin[i] = PETSC_NULL; 1534 1535 goto done; 1536 } 1537 } 1538 xs = dd->xs; 1539 ys = dd->ys; 1540 zs = dd->zs; 1541 xm = dd->xe-dd->xs; 1542 ym = dd->ye-dd->ys; 1543 zm = dd->ze-dd->zs; 1544 } 1545 1546 switch (dd->dim) { 1547 case 1: { 1548 void *ptr; 1549 itdof = xm; 1550 1551 ierr = PetscMalloc(xm*bs1*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1552 1553 ptr = (void*)(iarray_start - xs*bs1*sizeof(PetscScalar)); 1554 *iptr = (void*)ptr; 1555 break;} 1556 case 2: { 1557 void **ptr; 1558 itdof = xm*ym; 1559 1560 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*bs1*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1561 1562 ptr = (void**)(iarray_start + xm*ym*bs1*sizeof(PetscScalar) - ys*sizeof(void*)); 1563 for(j=ys;j<ys+ym;j++) { 1564 ptr[j] = iarray_start + bs1*sizeof(PetscScalar)*(xm*(j-ys) - xs); 1565 } 1566 *iptr = (void*)ptr; 1567 break;} 1568 case 3: { 1569 void ***ptr,**bptr; 1570 itdof = xm*ym*zm; 1571 1572 ierr = PetscMalloc((zm+1)*sizeof(void **)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*bs1*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1573 1574 ptr = (void***)(iarray_start + xm*ym*zm*2*sizeof(PetscScalar) - zs*sizeof(void*)); 1575 bptr = (void**)(iarray_start + xm*ym*zm*2*sizeof(PetscScalar) + zm*sizeof(void**)); 1576 for(i=zs;i<zs+zm;i++) { 1577 ptr[i] = bptr + ((i-zs)*ym* - ys)*sizeof(void*); 1578 } 1579 for (i=zs; i<zs+zm; i++) { 1580 for (j=ys; j<ys+ym; j++) { 1581 ptr[i][j] = iarray_start + bs1*sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs); 1582 } 1583 } 1584 1585 *iptr = (void*)ptr; 1586 break;} 1587 default: 1588 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 1589 } 1590 1591 done: 1592 /* add arrays to the checked out list */ 1593 if (ghosted) { 1594 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1595 if (!dd->admfarrayghostedout[i]) { 1596 dd->admfarrayghostedout[i] = *iptr ; 1597 dd->admfstartghostedout[i] = iarray_start; 1598 dd->ghostedtdof = itdof; 1599 break; 1600 } 1601 } 1602 } else { 1603 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1604 if (!dd->admfarrayout[i]) { 1605 dd->admfarrayout[i] = *iptr ; 1606 dd->admfstartout[i] = iarray_start; 1607 dd->tdof = itdof; 1608 break; 1609 } 1610 } 1611 } 1612 if (i == DMDA_MAX_AD_ARRAYS+1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Too many DMDA ADIC arrays obtained"); 1613 if (tdof) *tdof = itdof; 1614 if (array_start) *(void**)array_start = iarray_start; 1615 PetscFunctionReturn(0); 1616 } 1617 1618 #undef __FUNCT__ 1619 #define __FUNCT__ "DMDARestoreAdicMFArray" 1620 /*@C 1621 DMDARestoreAdicMFArray - Restores an array of derivative types for a DMDA. 1622 1623 Input Parameter: 1624 + da - information about my local patch 1625 - ghosted - do you want arrays for the ghosted or nonghosted patch? 1626 1627 Output Parameters: 1628 + ptr - array data structure to be passed to ad_FormFunctionLocal() 1629 . array_start - actual start of 1d array of all values that adiC can access directly 1630 - tdof - total number of degrees of freedom represented in array_start 1631 1632 Level: advanced 1633 1634 .seealso: DMDAGetAdicArray() 1635 1636 @*/ 1637 PetscErrorCode DMDARestoreAdicMFArray(DM da,PetscBool ghosted,void *vptr,void *array_start,PetscInt *tdof) 1638 { 1639 PetscInt i; 1640 void **iptr = (void**)vptr,*iarray_start = 0; 1641 DM_DA *dd = (DM_DA*)da->data; 1642 1643 PetscFunctionBegin; 1644 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1645 if (ghosted) { 1646 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1647 if (dd->admfarrayghostedout[i] == *iptr) { 1648 iarray_start = dd->admfstartghostedout[i]; 1649 dd->admfarrayghostedout[i] = PETSC_NULL; 1650 dd->admfstartghostedout[i] = PETSC_NULL; 1651 break; 1652 } 1653 } 1654 if (!iarray_start) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Could not find array in checkout list"); 1655 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1656 if (!dd->admfarrayghostedin[i]){ 1657 dd->admfarrayghostedin[i] = *iptr; 1658 dd->admfstartghostedin[i] = iarray_start; 1659 break; 1660 } 1661 } 1662 } else { 1663 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1664 if (dd->admfarrayout[i] == *iptr) { 1665 iarray_start = dd->admfstartout[i]; 1666 dd->admfarrayout[i] = PETSC_NULL; 1667 dd->admfstartout[i] = PETSC_NULL; 1668 break; 1669 } 1670 } 1671 if (!iarray_start) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Could not find array in checkout list"); 1672 for (i=0; i<DMDA_MAX_AD_ARRAYS; i++) { 1673 if (!dd->admfarrayin[i]){ 1674 dd->admfarrayin[i] = *iptr; 1675 dd->admfstartin[i] = iarray_start; 1676 break; 1677 } 1678 } 1679 } 1680 PetscFunctionReturn(0); 1681 } 1682 1683 #undef __FUNCT__ 1684 #define __FUNCT__ "admf_DAGetArray" 1685 PetscErrorCode admf_DAGetArray(DM da,PetscBool ghosted,void *iptr) 1686 { 1687 PetscErrorCode ierr; 1688 PetscFunctionBegin; 1689 ierr = DMDAGetAdicMFArray(da,ghosted,iptr,0,0);CHKERRQ(ierr); 1690 PetscFunctionReturn(0); 1691 } 1692 1693 #undef __FUNCT__ 1694 #define __FUNCT__ "admf_DARestoreArray" 1695 PetscErrorCode admf_DARestoreArray(DM da,PetscBool ghosted,void *iptr) 1696 { 1697 PetscErrorCode ierr; 1698 PetscFunctionBegin; 1699 ierr = DMDARestoreAdicMFArray(da,ghosted,iptr,0,0);CHKERRQ(ierr); 1700 PetscFunctionReturn(0); 1701 } 1702 1703