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