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 = VecGetDM(vec, &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 if (da->defaultSection) { 59 ierr = DMCreateLocalVector_Section_Private(da,g);CHKERRQ(ierr); 60 } else { 61 ierr = VecCreate(PETSC_COMM_SELF,g);CHKERRQ(ierr); 62 ierr = VecSetSizes(*g,dd->nlocal,PETSC_DETERMINE);CHKERRQ(ierr); 63 ierr = VecSetBlockSize(*g,dd->w);CHKERRQ(ierr); 64 ierr = VecSetType(*g,da->vectype);CHKERRQ(ierr); 65 ierr = VecSetDM(*g, da);CHKERRQ(ierr); 66 #if defined(PETSC_HAVE_MATLAB_ENGINE) 67 if (dd->w == 1 && dd->dim == 2) { 68 ierr = PetscObjectComposeFunctionDynamic((PetscObject)*g,"PetscMatlabEnginePut_C","VecMatlabEnginePut_DA2d",VecMatlabEnginePut_DA2d);CHKERRQ(ierr); 69 } 70 #endif 71 } 72 PetscFunctionReturn(0); 73 } 74 75 #undef __FUNCT__ 76 #define __FUNCT__ "DMDAGetNumCells" 77 PetscErrorCode DMDAGetNumCells(DM dm, PetscInt *numCells) 78 { 79 DM_DA *da = (DM_DA *) dm->data; 80 const PetscInt dim = da->dim; 81 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 82 const PetscInt nC = (mx)*(dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1); 83 84 PetscFunctionBegin; 85 if (numCells) { 86 PetscValidIntPointer(numCells,2); 87 *numCells = nC; 88 } 89 PetscFunctionReturn(0); 90 } 91 92 #undef __FUNCT__ 93 #define __FUNCT__ "DMDAGetNumVertices" 94 PetscErrorCode DMDAGetNumVertices(DM dm, PetscInt *numVerticesX, PetscInt *numVerticesY, PetscInt *numVerticesZ, PetscInt *numVertices) 95 { 96 DM_DA *da = (DM_DA *) dm->data; 97 const PetscInt dim = da->dim; 98 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 99 const PetscInt nVx = mx+1; 100 const PetscInt nVy = dim > 1 ? (my+1) : 1; 101 const PetscInt nVz = dim > 2 ? (mz+1) : 1; 102 const PetscInt nV = nVx*nVy*nVz; 103 104 PetscFunctionBegin; 105 if (numVerticesX) { 106 PetscValidIntPointer(numVerticesX,2); 107 *numVerticesX = nVx; 108 } 109 if (numVerticesY) { 110 PetscValidIntPointer(numVerticesY,3); 111 *numVerticesY = nVy; 112 } 113 if (numVerticesZ) { 114 PetscValidIntPointer(numVerticesZ,4); 115 *numVerticesZ = nVz; 116 } 117 if (numVertices) { 118 PetscValidIntPointer(numVertices,5); 119 *numVertices = nV; 120 } 121 PetscFunctionReturn(0); 122 } 123 124 #undef __FUNCT__ 125 #define __FUNCT__ "DMDAGetNumFaces" 126 PetscErrorCode DMDAGetNumFaces(DM dm, PetscInt *numXFacesX, PetscInt *numXFaces, PetscInt *numYFacesY, PetscInt *numYFaces, PetscInt *numZFacesZ, PetscInt *numZFaces) 127 { 128 DM_DA *da = (DM_DA *) dm->data; 129 const PetscInt dim = da->dim; 130 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 131 const PetscInt nxF = (dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1); 132 const PetscInt nXF = (mx+1)*nxF; 133 const PetscInt nyF = mx*(dim > 2 ? mz : 1); 134 const PetscInt nYF = dim > 1 ? (my+1)*nyF : 0; 135 const PetscInt nzF = mx*(dim > 1 ? my : 0); 136 const PetscInt nZF = dim > 2 ? (mz+1)*nzF : 0; 137 138 PetscFunctionBegin; 139 if (numXFacesX) { 140 PetscValidIntPointer(numXFacesX,2); 141 *numXFacesX = nxF; 142 } 143 if (numXFaces) { 144 PetscValidIntPointer(numXFaces,3); 145 *numXFaces = nXF; 146 } 147 if (numYFacesY) { 148 PetscValidIntPointer(numYFacesY,4); 149 *numYFacesY = nyF; 150 } 151 if (numYFaces) { 152 PetscValidIntPointer(numYFaces,5); 153 *numYFaces = nYF; 154 } 155 if (numZFacesZ) { 156 PetscValidIntPointer(numZFacesZ,6); 157 *numZFacesZ = nzF; 158 } 159 if (numZFaces) { 160 PetscValidIntPointer(numZFaces,7); 161 *numZFaces = nZF; 162 } 163 PetscFunctionReturn(0); 164 } 165 166 #undef __FUNCT__ 167 #define __FUNCT__ "DMDAGetHeightStratum" 168 PetscErrorCode DMDAGetHeightStratum(DM dm, PetscInt height, PetscInt *pStart, PetscInt *pEnd) 169 { 170 DM_DA *da = (DM_DA *) dm->data; 171 const PetscInt dim = da->dim; 172 PetscInt nC, nV, nXF, nYF, nZF; 173 PetscErrorCode ierr; 174 175 PetscFunctionBegin; 176 if (pStart) {PetscValidIntPointer(pStart,3);} 177 if (pEnd) {PetscValidIntPointer(pEnd,4);} 178 ierr = DMDAGetNumCells(dm, &nC);CHKERRQ(ierr); 179 ierr = DMDAGetNumVertices(dm, PETSC_NULL, PETSC_NULL, PETSC_NULL, &nV);CHKERRQ(ierr); 180 ierr = DMDAGetNumFaces(dm, PETSC_NULL, &nXF, PETSC_NULL, &nYF, PETSC_NULL, &nZF);CHKERRQ(ierr); 181 if (height == 0) { 182 /* Cells */ 183 if (pStart) {*pStart = 0;} 184 if (pEnd) {*pEnd = nC;} 185 } else if (height == 1) { 186 /* Faces */ 187 if (pStart) {*pStart = nC+nV;} 188 if (pEnd) {*pEnd = nC+nV+nXF+nYF+nZF;} 189 } else if (height == dim) { 190 /* Vertices */ 191 if (pStart) {*pStart = nC;} 192 if (pEnd) {*pEnd = nC+nV;} 193 } else if (height < 0) { 194 /* All points */ 195 if (pStart) {*pStart = 0;} 196 if (pEnd) {*pEnd = nC+nV+nXF+nYF+nZF;} 197 } else SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "No points of height %d in the DA", height); 198 PetscFunctionReturn(0); 199 } 200 201 #undef __FUNCT__ 202 #define __FUNCT__ "DMDACreateSection" 203 /*@C 204 DMDACreateSection - Create a PetscSection inside the DMDA that describes data layout. This allows multiple fields with 205 different numbers of dofs on vertices, cells, and faces in each direction. 206 207 Input Parameters: 208 + dm- The DMDA 209 . numFields - The number of fields 210 . numComp - The number of components in each field, or PETSC_NULL for 1 211 . numVertexDof - The number of dofs per vertex for each field, or PETSC_NULL 212 . numFaceDof - The number of dofs per face for each field and direction, or PETSC_NULL 213 - numCellDof - The number of dofs per cell for each field, or PETSC_NULL 214 215 Level: developer 216 217 Note: 218 The default DMDA numbering is as follows: 219 220 - Cells: [0, nC) 221 - Vertices: [nC, nC+nV) 222 - X-Faces: [nC+nV, nC+nV+nXF) normal is +- x-dir 223 - Y-Faces: [nC+nV+nXF, nC+nV+nXF+nYF) normal is +- y-dir 224 - Z-Faces: [nC+nV+nXF+nYF, nC+nV+nXF+nYF+nZF) normal is +- z-dir 225 226 We interpret the default DMDA partition as a cell partition, and the data assignment as a cell assignment. 227 @*/ 228 PetscErrorCode DMDACreateSection(DM dm, PetscInt numComp[], PetscInt numVertexDof[], PetscInt numFaceDof[], PetscInt numCellDof[]) 229 { 230 DM_DA *da = (DM_DA *) dm->data; 231 const PetscInt dim = da->dim; 232 PetscInt numFields, numVertexTotDof = 0, numCellTotDof = 0, numFaceTotDof[3] = {0, 0, 0}; 233 PetscSF sf; 234 PetscMPIInt rank; 235 const PetscMPIInt *neighbors; 236 PetscInt *localPoints; 237 PetscSFNode *remotePoints; 238 PetscInt nleaves = 0, nleavesCheck = 0, nL = 0; 239 PetscInt nC, nVx, nVy, nVz, nV, nxF, nXF, nyF, nYF, nzF, nZF; 240 PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart, zfEnd; 241 PetscInt f, v, c, xf, yf, zf, xn, yn, zn; 242 PetscErrorCode ierr; 243 244 PetscFunctionBegin; 245 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 246 ierr = MPI_Comm_rank(((PetscObject) dm)->comm, &rank);CHKERRQ(ierr); 247 ierr = DMDAGetNumCells(dm, &nC);CHKERRQ(ierr); 248 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr); 249 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 250 ierr = DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);CHKERRQ(ierr); 251 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 252 ierr = DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 253 ierr = DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);CHKERRQ(ierr); 254 xfStart = vEnd; xfEnd = xfStart+nXF; 255 yfStart = xfEnd; yfEnd = yfStart+nYF; 256 zfStart = yfEnd; zfEnd = zfStart+nZF; 257 if (zfEnd != fEnd) SETERRQ2(((PetscObject) dm)->comm, PETSC_ERR_PLIB, "Invalid face end %d, should be %d", zfEnd, fEnd); 258 /* Create local section */ 259 ierr = DMDAGetInfo(dm, 0,0,0,0,0,0,0, &numFields, 0,0,0,0,0);CHKERRQ(ierr); 260 for (f = 0; f < numFields; ++f) { 261 if (numVertexDof) {numVertexTotDof += numVertexDof[f];} 262 if (numCellDof) {numCellTotDof += numCellDof[f];} 263 if (numFaceDof) {numFaceTotDof[0] += numFaceDof[f*dim+0]; 264 numFaceTotDof[1] += dim > 1 ? numFaceDof[f*dim+1] : 0; 265 numFaceTotDof[2] += dim > 2 ? numFaceDof[f*dim+2] : 0; 266 } 267 } 268 ierr = PetscSectionCreate(((PetscObject) dm)->comm, &dm->defaultSection);CHKERRQ(ierr); 269 if (numFields > 1) { 270 ierr = PetscSectionSetNumFields(dm->defaultSection, numFields);CHKERRQ(ierr); 271 for (f = 0; f < numFields; ++f) { 272 const char *name; 273 274 ierr = DMDAGetFieldName(dm, f, &name);CHKERRQ(ierr); 275 ierr = PetscSectionSetFieldName(dm->defaultSection, f, name);CHKERRQ(ierr); 276 if (numComp) { 277 ierr = PetscSectionSetFieldComponents(dm->defaultSection, f, numComp[f]);CHKERRQ(ierr); 278 } 279 } 280 } else { 281 numFields = 0; 282 } 283 ierr = PetscSectionSetChart(dm->defaultSection, pStart, pEnd);CHKERRQ(ierr); 284 if (numVertexDof) { 285 for (v = vStart; v < vEnd; ++v) { 286 for (f = 0; f < numFields; ++f) { 287 ierr = PetscSectionSetFieldDof(dm->defaultSection, v, f, numVertexDof[f]);CHKERRQ(ierr); 288 } 289 ierr = PetscSectionSetDof(dm->defaultSection, v, numVertexTotDof);CHKERRQ(ierr); 290 } 291 } 292 if (numFaceDof) { 293 for (xf = xfStart; xf < xfEnd; ++xf) { 294 for (f = 0; f < numFields; ++f) { 295 ierr = PetscSectionSetFieldDof(dm->defaultSection, xf, f, numFaceDof[f*dim+0]);CHKERRQ(ierr); 296 } 297 ierr = PetscSectionSetDof(dm->defaultSection, xf, numFaceTotDof[0]);CHKERRQ(ierr); 298 } 299 for (yf = yfStart; yf < yfEnd; ++yf) { 300 for (f = 0; f < numFields; ++f) { 301 ierr = PetscSectionSetFieldDof(dm->defaultSection, yf, f, numFaceDof[f*dim+1]);CHKERRQ(ierr); 302 } 303 ierr = PetscSectionSetDof(dm->defaultSection, yf, numFaceTotDof[1]);CHKERRQ(ierr); 304 } 305 for (zf = zfStart; zf < zfEnd; ++zf) { 306 for (f = 0; f < numFields; ++f) { 307 ierr = PetscSectionSetFieldDof(dm->defaultSection, zf, f, numFaceDof[f*dim+2]);CHKERRQ(ierr); 308 } 309 ierr = PetscSectionSetDof(dm->defaultSection, zf, numFaceTotDof[2]);CHKERRQ(ierr); 310 } 311 } 312 if (numCellDof) { 313 for (c = cStart; c < cEnd; ++c) { 314 for (f = 0; f < numFields; ++f) { 315 ierr = PetscSectionSetFieldDof(dm->defaultSection, c, f, numCellDof[f]);CHKERRQ(ierr); 316 } 317 ierr = PetscSectionSetDof(dm->defaultSection, c, numCellTotDof);CHKERRQ(ierr); 318 } 319 } 320 ierr = PetscSectionSetUp(dm->defaultSection);CHKERRQ(ierr); 321 /* Create mesh point SF */ 322 ierr = DMDAGetNeighbors(dm, &neighbors);CHKERRQ(ierr); 323 for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) { 324 for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) { 325 for (xn = 0; xn < 3; ++xn) { 326 const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0; 327 const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn]; 328 329 if (neighbor >= 0 && neighbor < rank) { 330 nleaves += (!xp ? nVx : 1) * (!yp ? nVy : 1) * (!zp ? nVz : 1); /* vertices */ 331 if (xp && !yp && !zp) { 332 nleaves += nxF; /* x faces */ 333 } else if (yp && !zp && !xp) { 334 nleaves += nyF; /* y faces */ 335 } else if (zp && !xp && !yp) { 336 nleaves += nzF; /* z faces */ 337 } 338 } 339 } 340 } 341 } 342 ierr = PetscMalloc2(nleaves,PetscInt,&localPoints,nleaves,PetscSFNode,&remotePoints);CHKERRQ(ierr); 343 for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) { 344 for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) { 345 for (xn = 0; xn < 3; ++xn) { 346 const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0; 347 const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn]; 348 PetscInt xv, yv, zv; 349 350 if (neighbor >= 0 && neighbor < rank) { 351 if (xp < 0) { /* left */ 352 if (yp < 0) { /* bottom */ 353 if (zp < 0) { /* back */ 354 const PetscInt localVertex = ( 0*nVy + 0)*nVx + 0 + nC; 355 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 356 nleavesCheck += 1; /* left bottom back vertex */ 357 358 localPoints[nL] = localVertex; 359 remotePoints[nL].rank = neighbor; 360 remotePoints[nL].index = remoteVertex; 361 ++nL; 362 } else if (zp > 0) { /* front */ 363 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; 364 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 365 nleavesCheck += 1; /* left bottom front vertex */ 366 367 localPoints[nL] = localVertex; 368 remotePoints[nL].rank = neighbor; 369 remotePoints[nL].index = remoteVertex; 370 ++nL; 371 } else { 372 nleavesCheck += nVz; /* left bottom vertices */ 373 for (zv = 0; zv < nVz; ++zv, ++nL) { 374 const PetscInt localVertex = (zv*nVy + 0)*nVx + 0 + nC; 375 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 376 377 localPoints[nL] = localVertex; 378 remotePoints[nL].rank = neighbor; 379 remotePoints[nL].index = remoteVertex; 380 } 381 } 382 } else if (yp > 0) { /* top */ 383 if (zp < 0) { /* back */ 384 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; 385 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 386 nleavesCheck += 1; /* left top back vertex */ 387 388 localPoints[nL] = localVertex; 389 remotePoints[nL].rank = neighbor; 390 remotePoints[nL].index = remoteVertex; 391 ++nL; 392 } else if (zp > 0) { /* front */ 393 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; 394 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 395 nleavesCheck += 1; /* left top front vertex */ 396 397 localPoints[nL] = localVertex; 398 remotePoints[nL].rank = neighbor; 399 remotePoints[nL].index = remoteVertex; 400 ++nL; 401 } else { 402 nleavesCheck += nVz; /* left top vertices */ 403 for (zv = 0; zv < nVz; ++zv, ++nL) { 404 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; 405 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 406 407 localPoints[nL] = localVertex; 408 remotePoints[nL].rank = neighbor; 409 remotePoints[nL].index = remoteVertex; 410 } 411 } 412 } else { 413 if (zp < 0) { /* back */ 414 nleavesCheck += nVy; /* left back vertices */ 415 for (yv = 0; yv < nVy; ++yv, ++nL) { 416 const PetscInt localVertex = ( 0*nVy + yv)*nVx + 0 + nC; 417 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 418 419 localPoints[nL] = localVertex; 420 remotePoints[nL].rank = neighbor; 421 remotePoints[nL].index = remoteVertex; 422 } 423 } else if (zp > 0) { /* front */ 424 nleavesCheck += nVy; /* left front vertices */ 425 for (yv = 0; yv < nVy; ++yv, ++nL) { 426 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; 427 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 428 429 localPoints[nL] = localVertex; 430 remotePoints[nL].rank = neighbor; 431 remotePoints[nL].index = remoteVertex; 432 } 433 } else { 434 nleavesCheck += nVy*nVz; /* left vertices */ 435 for (zv = 0; zv < nVz; ++zv) { 436 for (yv = 0; yv < nVy; ++yv, ++nL) { 437 const PetscInt localVertex = (zv*nVy + yv)*nVx + 0 + nC; 438 const PetscInt remoteVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 439 440 localPoints[nL] = localVertex; 441 remotePoints[nL].rank = neighbor; 442 remotePoints[nL].index = remoteVertex; 443 } 444 } 445 nleavesCheck += nxF; /* left faces */ 446 for (xf = 0; xf < nxF; ++xf, ++nL) { 447 /* THIS IS WRONG */ 448 const PetscInt localFace = 0 + nC+nV; 449 const PetscInt remoteFace = 0 + nC+nV; 450 451 localPoints[nL] = localFace; 452 remotePoints[nL].rank = neighbor; 453 remotePoints[nL].index = remoteFace; 454 } 455 } 456 } 457 } else if (xp > 0) { /* right */ 458 if (yp < 0) { /* bottom */ 459 if (zp < 0) { /* back */ 460 const PetscInt localVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; 461 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 462 nleavesCheck += 1; /* right bottom back vertex */ 463 464 localPoints[nL] = localVertex; 465 remotePoints[nL].rank = neighbor; 466 remotePoints[nL].index = remoteVertex; 467 ++nL; 468 } else if (zp > 0) { /* front */ 469 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; 470 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 471 nleavesCheck += 1; /* right bottom front vertex */ 472 473 localPoints[nL] = localVertex; 474 remotePoints[nL].rank = neighbor; 475 remotePoints[nL].index = remoteVertex; 476 ++nL; 477 } else { 478 nleavesCheck += nVz; /* right bottom vertices */ 479 for (zv = 0; zv < nVz; ++zv, ++nL) { 480 const PetscInt localVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; 481 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 482 483 localPoints[nL] = localVertex; 484 remotePoints[nL].rank = neighbor; 485 remotePoints[nL].index = remoteVertex; 486 } 487 } 488 } else if (yp > 0) { /* top */ 489 if (zp < 0) { /* back */ 490 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; 491 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 492 nleavesCheck += 1; /* right top back vertex */ 493 494 localPoints[nL] = localVertex; 495 remotePoints[nL].rank = neighbor; 496 remotePoints[nL].index = remoteVertex; 497 ++nL; 498 } else if (zp > 0) { /* front */ 499 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; 500 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 501 nleavesCheck += 1; /* right top front vertex */ 502 503 localPoints[nL] = localVertex; 504 remotePoints[nL].rank = neighbor; 505 remotePoints[nL].index = remoteVertex; 506 ++nL; 507 } else { 508 nleavesCheck += nVz; /* right top vertices */ 509 for (zv = 0; zv < nVz; ++zv, ++nL) { 510 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; 511 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 512 513 localPoints[nL] = localVertex; 514 remotePoints[nL].rank = neighbor; 515 remotePoints[nL].index = remoteVertex; 516 } 517 } 518 } else { 519 if (zp < 0) { /* back */ 520 nleavesCheck += nVy; /* right back vertices */ 521 for (yv = 0; yv < nVy; ++yv, ++nL) { 522 const PetscInt localVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; 523 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 524 525 localPoints[nL] = localVertex; 526 remotePoints[nL].rank = neighbor; 527 remotePoints[nL].index = remoteVertex; 528 } 529 } else if (zp > 0) { /* front */ 530 nleavesCheck += nVy; /* right front vertices */ 531 for (yv = 0; yv < nVy; ++yv, ++nL) { 532 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; 533 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 534 535 localPoints[nL] = localVertex; 536 remotePoints[nL].rank = neighbor; 537 remotePoints[nL].index = remoteVertex; 538 } 539 } else { 540 nleavesCheck += nVy*nVz; /* right vertices */ 541 for (zv = 0; zv < nVz; ++zv) { 542 for (yv = 0; yv < nVy; ++yv, ++nL) { 543 const PetscInt localVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; 544 const PetscInt remoteVertex = (zv*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 545 546 localPoints[nL] = localVertex; 547 remotePoints[nL].rank = neighbor; 548 remotePoints[nL].index = remoteVertex; 549 } 550 } 551 nleavesCheck += nxF; /* right faces */ 552 for (xf = 0; xf < nxF; ++xf, ++nL) { 553 /* THIS IS WRONG */ 554 const PetscInt localFace = 0 + nC+nV; 555 const PetscInt remoteFace = 0 + nC+nV; 556 557 localPoints[nL] = localFace; 558 remotePoints[nL].rank = neighbor; 559 remotePoints[nL].index = remoteFace; 560 } 561 } 562 } 563 } else { 564 if (yp < 0) { /* bottom */ 565 if (zp < 0) { /* back */ 566 nleavesCheck += nVx; /* bottom back vertices */ 567 for (xv = 0; xv < nVx; ++xv, ++nL) { 568 const PetscInt localVertex = ( 0*nVy + 0)*nVx + xv + nC; 569 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 570 571 localPoints[nL] = localVertex; 572 remotePoints[nL].rank = neighbor; 573 remotePoints[nL].index = remoteVertex; 574 } 575 } else if (zp > 0) { /* front */ 576 nleavesCheck += nVx; /* bottom front vertices */ 577 for (xv = 0; xv < nVx; ++xv, ++nL) { 578 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; 579 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 580 581 localPoints[nL] = localVertex; 582 remotePoints[nL].rank = neighbor; 583 remotePoints[nL].index = remoteVertex; 584 } 585 } else { 586 nleavesCheck += nVx*nVz; /* bottom vertices */ 587 for (zv = 0; zv < nVz; ++zv) { 588 for (xv = 0; xv < nVx; ++xv, ++nL) { 589 const PetscInt localVertex = (zv*nVy + 0)*nVx + xv + nC; 590 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 591 592 localPoints[nL] = localVertex; 593 remotePoints[nL].rank = neighbor; 594 remotePoints[nL].index = remoteVertex; 595 } 596 } 597 nleavesCheck += nyF; /* bottom faces */ 598 for (yf = 0; yf < nyF; ++yf, ++nL) { 599 /* THIS IS WRONG */ 600 const PetscInt localFace = 0 + nC+nV; 601 const PetscInt remoteFace = 0 + nC+nV; 602 603 localPoints[nL] = localFace; 604 remotePoints[nL].rank = neighbor; 605 remotePoints[nL].index = remoteFace; 606 } 607 } 608 } else if (yp > 0) { /* top */ 609 if (zp < 0) { /* back */ 610 nleavesCheck += nVx; /* top back vertices */ 611 for (xv = 0; xv < nVx; ++xv, ++nL) { 612 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; 613 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 614 615 localPoints[nL] = localVertex; 616 remotePoints[nL].rank = neighbor; 617 remotePoints[nL].index = remoteVertex; 618 } 619 } else if (zp > 0) { /* front */ 620 nleavesCheck += nVx; /* top front vertices */ 621 for (xv = 0; xv < nVx; ++xv, ++nL) { 622 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; 623 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 624 625 localPoints[nL] = localVertex; 626 remotePoints[nL].rank = neighbor; 627 remotePoints[nL].index = remoteVertex; 628 } 629 } else { 630 nleavesCheck += nVx*nVz; /* top vertices */ 631 for (zv = 0; zv < nVz; ++zv) { 632 for (xv = 0; xv < nVx; ++xv, ++nL) { 633 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + xv + nC; 634 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 635 636 localPoints[nL] = localVertex; 637 remotePoints[nL].rank = neighbor; 638 remotePoints[nL].index = remoteVertex; 639 } 640 } 641 nleavesCheck += nyF; /* top faces */ 642 for (yf = 0; yf < nyF; ++yf, ++nL) { 643 /* THIS IS WRONG */ 644 const PetscInt localFace = 0 + nC+nV; 645 const PetscInt remoteFace = 0 + nC+nV; 646 647 localPoints[nL] = localFace; 648 remotePoints[nL].rank = neighbor; 649 remotePoints[nL].index = remoteFace; 650 } 651 } 652 } else { 653 if (zp < 0) { /* back */ 654 nleavesCheck += nVx*nVy; /* back vertices */ 655 for (yv = 0; yv < nVy; ++yv) { 656 for (xv = 0; xv < nVx; ++xv, ++nL) { 657 const PetscInt localVertex = ( 0*nVy + yv)*nVx + xv + nC; 658 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 659 660 localPoints[nL] = localVertex; 661 remotePoints[nL].rank = neighbor; 662 remotePoints[nL].index = remoteVertex; 663 } 664 } 665 nleavesCheck += nzF; /* back faces */ 666 for (zf = 0; zf < nzF; ++zf, ++nL) { 667 /* THIS IS WRONG */ 668 const PetscInt localFace = 0 + nC+nV; 669 const PetscInt remoteFace = 0 + nC+nV; 670 671 localPoints[nL] = localFace; 672 remotePoints[nL].rank = neighbor; 673 remotePoints[nL].index = remoteFace; 674 } 675 } else if (zp > 0) { /* front */ 676 nleavesCheck += nVx*nVy; /* front vertices */ 677 for (yv = 0; yv < nVy; ++yv) { 678 for (xv = 0; xv < nVx; ++xv, ++nL) { 679 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; 680 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 681 682 localPoints[nL] = localVertex; 683 remotePoints[nL].rank = neighbor; 684 remotePoints[nL].index = remoteVertex; 685 } 686 } 687 nleavesCheck += nzF; /* front faces */ 688 for (zf = 0; zf < nzF; ++zf, ++nL) { 689 /* THIS IS WRONG */ 690 const PetscInt localFace = 0 + nC+nV; 691 const PetscInt remoteFace = 0 + nC+nV; 692 693 localPoints[nL] = localFace; 694 remotePoints[nL].rank = neighbor; 695 remotePoints[nL].index = remoteFace; 696 } 697 } else { 698 /* Nothing is shared from the interior */ 699 } 700 } 701 } 702 } 703 } 704 } 705 } 706 /* TODO: Remove duplication in leaf determination */ 707 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); 708 ierr = PetscSFCreate(((PetscObject) dm)->comm, &sf);CHKERRQ(ierr); 709 ierr = PetscSFSetGraph(sf, pEnd, nleaves, localPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 710 /* Create global section */ 711 ierr = PetscSectionCreateGlobalSection(dm->defaultSection, sf, PETSC_FALSE, &dm->defaultGlobalSection);CHKERRQ(ierr); 712 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 713 /* Create default SF */ 714 ierr = DMCreateDefaultSF(dm, dm->defaultSection, dm->defaultGlobalSection);CHKERRQ(ierr); 715 PetscFunctionReturn(0); 716 } 717 718 /* ------------------------------------------------------------------- */ 719 720 #undef __FUNCT__ 721 #define __FUNCT__ "DMDAGetArray" 722 /*@C 723 DMDAGetArray - Gets a work array for a DMDA 724 725 Input Parameter: 726 + da - information about my local patch 727 - ghosted - do you want arrays for the ghosted or nonghosted patch 728 729 Output Parameters: 730 . vptr - array data structured 731 732 Note: The vector values are NOT initialized and may have garbage in them, so you may need 733 to zero them. 734 735 Level: advanced 736 737 .seealso: DMDARestoreArray() 738 739 @*/ 740 PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr) 741 { 742 PetscErrorCode ierr; 743 PetscInt j,i,xs,ys,xm,ym,zs,zm; 744 char *iarray_start; 745 void **iptr = (void**)vptr; 746 DM_DA *dd = (DM_DA*)da->data; 747 748 PetscFunctionBegin; 749 PetscValidHeaderSpecific(da,DM_CLASSID,1); 750 if (ghosted) { 751 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 752 if (dd->arrayghostedin[i]) { 753 *iptr = dd->arrayghostedin[i]; 754 iarray_start = (char*)dd->startghostedin[i]; 755 dd->arrayghostedin[i] = PETSC_NULL; 756 dd->startghostedin[i] = PETSC_NULL; 757 758 goto done; 759 } 760 } 761 xs = dd->Xs; 762 ys = dd->Ys; 763 zs = dd->Zs; 764 xm = dd->Xe-dd->Xs; 765 ym = dd->Ye-dd->Ys; 766 zm = dd->Ze-dd->Zs; 767 } else { 768 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 769 if (dd->arrayin[i]) { 770 *iptr = dd->arrayin[i]; 771 iarray_start = (char*)dd->startin[i]; 772 dd->arrayin[i] = PETSC_NULL; 773 dd->startin[i] = PETSC_NULL; 774 775 goto done; 776 } 777 } 778 xs = dd->xs; 779 ys = dd->ys; 780 zs = dd->zs; 781 xm = dd->xe-dd->xs; 782 ym = dd->ye-dd->ys; 783 zm = dd->ze-dd->zs; 784 } 785 786 switch (dd->dim) { 787 case 1: { 788 void *ptr; 789 790 ierr = PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 791 792 ptr = (void*)(iarray_start - xs*sizeof(PetscScalar)); 793 *iptr = (void*)ptr; 794 break;} 795 case 2: { 796 void **ptr; 797 798 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 799 800 ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*)); 801 for (j=ys;j<ys+ym;j++) { 802 ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs); 803 } 804 *iptr = (void*)ptr; 805 break;} 806 case 3: { 807 void ***ptr,**bptr; 808 809 ierr = PetscMalloc((zm+1)*sizeof(void **)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 810 811 ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*)); 812 bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**)); 813 for (i=zs;i<zs+zm;i++) { 814 ptr[i] = bptr + ((i-zs)*ym - ys); 815 } 816 for (i=zs; i<zs+zm; i++) { 817 for (j=ys; j<ys+ym; j++) { 818 ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs); 819 } 820 } 821 822 *iptr = (void*)ptr; 823 break;} 824 default: 825 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 826 } 827 828 done: 829 /* add arrays to the checked out list */ 830 if (ghosted) { 831 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 832 if (!dd->arrayghostedout[i]) { 833 dd->arrayghostedout[i] = *iptr ; 834 dd->startghostedout[i] = iarray_start; 835 break; 836 } 837 } 838 } else { 839 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 840 if (!dd->arrayout[i]) { 841 dd->arrayout[i] = *iptr ; 842 dd->startout[i] = iarray_start; 843 break; 844 } 845 } 846 } 847 PetscFunctionReturn(0); 848 } 849 850 #undef __FUNCT__ 851 #define __FUNCT__ "DMDARestoreArray" 852 /*@C 853 DMDARestoreArray - Restores an array of derivative types for a DMDA 854 855 Input Parameter: 856 + da - information about my local patch 857 . ghosted - do you want arrays for the ghosted or nonghosted patch 858 - vptr - array data structured to be passed to ad_FormFunctionLocal() 859 860 Level: advanced 861 862 .seealso: DMDAGetArray() 863 864 @*/ 865 PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr) 866 { 867 PetscInt i; 868 void **iptr = (void**)vptr,*iarray_start = 0; 869 DM_DA *dd = (DM_DA*)da->data; 870 871 PetscFunctionBegin; 872 PetscValidHeaderSpecific(da,DM_CLASSID,1); 873 if (ghosted) { 874 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 875 if (dd->arrayghostedout[i] == *iptr) { 876 iarray_start = dd->startghostedout[i]; 877 dd->arrayghostedout[i] = PETSC_NULL; 878 dd->startghostedout[i] = PETSC_NULL; 879 break; 880 } 881 } 882 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 883 if (!dd->arrayghostedin[i]) { 884 dd->arrayghostedin[i] = *iptr; 885 dd->startghostedin[i] = iarray_start; 886 break; 887 } 888 } 889 } else { 890 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 891 if (dd->arrayout[i] == *iptr) { 892 iarray_start = dd->startout[i]; 893 dd->arrayout[i] = PETSC_NULL; 894 dd->startout[i] = PETSC_NULL; 895 break; 896 } 897 } 898 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 899 if (!dd->arrayin[i]) { 900 dd->arrayin[i] = *iptr; 901 dd->startin[i] = iarray_start; 902 break; 903 } 904 } 905 } 906 PetscFunctionReturn(0); 907 } 908 909