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 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 = VecSetDM(*g, 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 const PetscInt localVertex = ( 0*nVy + 0)*nVx + 0 + nC; 352 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 353 nleavesCheck += 1; /* left bottom back vertex */ 354 355 localPoints[nL] = localVertex; 356 remotePoints[nL].rank = neighbor; 357 remotePoints[nL].index = remoteVertex; 358 ++nL; 359 } else if (zp > 0) { /* front */ 360 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; 361 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 362 nleavesCheck += 1; /* left bottom front vertex */ 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 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; 382 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 383 nleavesCheck += 1; /* left top back vertex */ 384 385 localPoints[nL] = localVertex; 386 remotePoints[nL].rank = neighbor; 387 remotePoints[nL].index = remoteVertex; 388 ++nL; 389 } else if (zp > 0) { /* front */ 390 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; 391 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 392 nleavesCheck += 1; /* left top front vertex */ 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 const PetscInt localVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; 458 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 459 nleavesCheck += 1; /* right bottom back vertex */ 460 461 localPoints[nL] = localVertex; 462 remotePoints[nL].rank = neighbor; 463 remotePoints[nL].index = remoteVertex; 464 ++nL; 465 } else if (zp > 0) { /* front */ 466 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; 467 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 468 nleavesCheck += 1; /* right bottom front vertex */ 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 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; 488 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 489 nleavesCheck += 1; /* right top back vertex */ 490 491 localPoints[nL] = localVertex; 492 remotePoints[nL].rank = neighbor; 493 remotePoints[nL].index = remoteVertex; 494 ++nL; 495 } else if (zp > 0) { /* front */ 496 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; 497 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 498 nleavesCheck += 1; /* right top front vertex */ 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, PETSC_FALSE, &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 717 #undef __FUNCT__ 718 #define __FUNCT__ "DMDAGetArray" 719 /*@C 720 DMDAGetArray - Gets a work array for a DMDA 721 722 Input Parameter: 723 + da - information about my local patch 724 - ghosted - do you want arrays for the ghosted or nonghosted patch 725 726 Output Parameters: 727 . vptr - array data structured 728 729 Note: The vector values are NOT initialized and may have garbage in them, so you may need 730 to zero them. 731 732 Level: advanced 733 734 .seealso: DMDARestoreArray() 735 736 @*/ 737 PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr) 738 { 739 PetscErrorCode ierr; 740 PetscInt j,i,xs,ys,xm,ym,zs,zm; 741 char *iarray_start; 742 void **iptr = (void**)vptr; 743 DM_DA *dd = (DM_DA*)da->data; 744 745 PetscFunctionBegin; 746 PetscValidHeaderSpecific(da,DM_CLASSID,1); 747 if (ghosted) { 748 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 749 if (dd->arrayghostedin[i]) { 750 *iptr = dd->arrayghostedin[i]; 751 iarray_start = (char*)dd->startghostedin[i]; 752 dd->arrayghostedin[i] = PETSC_NULL; 753 dd->startghostedin[i] = PETSC_NULL; 754 755 goto done; 756 } 757 } 758 xs = dd->Xs; 759 ys = dd->Ys; 760 zs = dd->Zs; 761 xm = dd->Xe-dd->Xs; 762 ym = dd->Ye-dd->Ys; 763 zm = dd->Ze-dd->Zs; 764 } else { 765 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 766 if (dd->arrayin[i]) { 767 *iptr = dd->arrayin[i]; 768 iarray_start = (char*)dd->startin[i]; 769 dd->arrayin[i] = PETSC_NULL; 770 dd->startin[i] = PETSC_NULL; 771 772 goto done; 773 } 774 } 775 xs = dd->xs; 776 ys = dd->ys; 777 zs = dd->zs; 778 xm = dd->xe-dd->xs; 779 ym = dd->ye-dd->ys; 780 zm = dd->ze-dd->zs; 781 } 782 783 switch (dd->dim) { 784 case 1: { 785 void *ptr; 786 787 ierr = PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 788 789 ptr = (void*)(iarray_start - xs*sizeof(PetscScalar)); 790 *iptr = (void*)ptr; 791 break;} 792 case 2: { 793 void **ptr; 794 795 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 796 797 ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*)); 798 for (j=ys;j<ys+ym;j++) { 799 ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs); 800 } 801 *iptr = (void*)ptr; 802 break;} 803 case 3: { 804 void ***ptr,**bptr; 805 806 ierr = PetscMalloc((zm+1)*sizeof(void **)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 807 808 ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*)); 809 bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**)); 810 for (i=zs;i<zs+zm;i++) { 811 ptr[i] = bptr + ((i-zs)*ym - ys); 812 } 813 for (i=zs; i<zs+zm; i++) { 814 for (j=ys; j<ys+ym; j++) { 815 ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs); 816 } 817 } 818 819 *iptr = (void*)ptr; 820 break;} 821 default: 822 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 823 } 824 825 done: 826 /* add arrays to the checked out list */ 827 if (ghosted) { 828 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 829 if (!dd->arrayghostedout[i]) { 830 dd->arrayghostedout[i] = *iptr ; 831 dd->startghostedout[i] = iarray_start; 832 break; 833 } 834 } 835 } else { 836 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 837 if (!dd->arrayout[i]) { 838 dd->arrayout[i] = *iptr ; 839 dd->startout[i] = iarray_start; 840 break; 841 } 842 } 843 } 844 PetscFunctionReturn(0); 845 } 846 847 #undef __FUNCT__ 848 #define __FUNCT__ "DMDARestoreArray" 849 /*@C 850 DMDARestoreArray - Restores an array of derivative types for a DMDA 851 852 Input Parameter: 853 + da - information about my local patch 854 . ghosted - do you want arrays for the ghosted or nonghosted patch 855 - vptr - array data structured to be passed to ad_FormFunctionLocal() 856 857 Level: advanced 858 859 .seealso: DMDAGetArray() 860 861 @*/ 862 PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr) 863 { 864 PetscInt i; 865 void **iptr = (void**)vptr,*iarray_start = 0; 866 DM_DA *dd = (DM_DA*)da->data; 867 868 PetscFunctionBegin; 869 PetscValidHeaderSpecific(da,DM_CLASSID,1); 870 if (ghosted) { 871 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 872 if (dd->arrayghostedout[i] == *iptr) { 873 iarray_start = dd->startghostedout[i]; 874 dd->arrayghostedout[i] = PETSC_NULL; 875 dd->startghostedout[i] = PETSC_NULL; 876 break; 877 } 878 } 879 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 880 if (!dd->arrayghostedin[i]){ 881 dd->arrayghostedin[i] = *iptr; 882 dd->startghostedin[i] = iarray_start; 883 break; 884 } 885 } 886 } else { 887 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 888 if (dd->arrayout[i] == *iptr) { 889 iarray_start = dd->startout[i]; 890 dd->arrayout[i] = PETSC_NULL; 891 dd->startout[i] = PETSC_NULL; 892 break; 893 } 894 } 895 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 896 if (!dd->arrayin[i]){ 897 dd->arrayin[i] = *iptr; 898 dd->startin[i] = iarray_start; 899 break; 900 } 901 } 902 } 903 PetscFunctionReturn(0); 904 } 905 906