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