1 2 /* 3 Code for manipulating distributed regular arrays in parallel. 4 */ 5 6 #include <petsc-private/dmdaimpl.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 #undef __FUNCT__ 15 #define __FUNCT__ "VecMatlabEnginePut_DA2d" 16 static PetscErrorCode VecMatlabEnginePut_DA2d(PetscObject obj,void *mengine) 17 { 18 PetscErrorCode ierr; 19 PetscInt n,m; 20 Vec vec = (Vec)obj; 21 PetscScalar *array; 22 mxArray *mat; 23 DM da; 24 25 PetscFunctionBegin; 26 ierr = VecGetDM(vec, &da);CHKERRQ(ierr); 27 if (!da) SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_ARG_WRONGSTATE,"Vector not associated with a DMDA"); 28 ierr = DMDAGetGhostCorners(da,0,0,0,&m,&n,0);CHKERRQ(ierr); 29 30 ierr = VecGetArray(vec,&array);CHKERRQ(ierr); 31 #if !defined(PETSC_USE_COMPLEX) 32 mat = mxCreateDoubleMatrix(m,n,mxREAL); 33 #else 34 mat = mxCreateDoubleMatrix(m,n,mxCOMPLEX); 35 #endif 36 ierr = PetscMemcpy(mxGetPr(mat),array,n*m*sizeof(PetscScalar));CHKERRQ(ierr); 37 ierr = PetscObjectName(obj);CHKERRQ(ierr); 38 engPutVariable((Engine*)mengine,obj->name,mat); 39 40 ierr = VecRestoreArray(vec,&array);CHKERRQ(ierr); 41 PetscFunctionReturn(0); 42 } 43 #endif 44 45 46 #undef __FUNCT__ 47 #define __FUNCT__ "DMCreateLocalVector_DA" 48 PetscErrorCode DMCreateLocalVector_DA(DM da,Vec *g) 49 { 50 PetscErrorCode ierr; 51 DM_DA *dd = (DM_DA*)da->data; 52 53 PetscFunctionBegin; 54 PetscValidHeaderSpecific(da,DM_CLASSID,1); 55 PetscValidPointer(g,2); 56 if (da->defaultSection) { 57 ierr = DMCreateLocalVector_Section_Private(da,g);CHKERRQ(ierr); 58 } else { 59 ierr = VecCreate(PETSC_COMM_SELF,g);CHKERRQ(ierr); 60 ierr = VecSetSizes(*g,dd->nlocal,PETSC_DETERMINE);CHKERRQ(ierr); 61 ierr = VecSetBlockSize(*g,dd->w);CHKERRQ(ierr); 62 ierr = VecSetType(*g,da->vectype);CHKERRQ(ierr); 63 ierr = VecSetDM(*g, da);CHKERRQ(ierr); 64 #if defined(PETSC_HAVE_MATLAB_ENGINE) 65 if (dd->w == 1 && dd->dim == 2) { 66 ierr = PetscObjectComposeFunction((PetscObject)*g,"PetscMatlabEnginePut_C",VecMatlabEnginePut_DA2d);CHKERRQ(ierr); 67 } 68 #endif 69 } 70 PetscFunctionReturn(0); 71 } 72 73 #undef __FUNCT__ 74 #define __FUNCT__ "DMDAGetNumCells" 75 PetscErrorCode DMDAGetNumCells(DM dm, PetscInt *numCells) 76 { 77 DM_DA *da = (DM_DA*) dm->data; 78 const PetscInt dim = da->dim; 79 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 80 const PetscInt nC = (mx)*(dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1); 81 82 PetscFunctionBegin; 83 if (numCells) { 84 PetscValidIntPointer(numCells,2); 85 *numCells = nC; 86 } 87 PetscFunctionReturn(0); 88 } 89 90 #undef __FUNCT__ 91 #define __FUNCT__ "DMDAGetNumVertices" 92 PetscErrorCode DMDAGetNumVertices(DM dm, PetscInt *numVerticesX, PetscInt *numVerticesY, PetscInt *numVerticesZ, PetscInt *numVertices) 93 { 94 DM_DA *da = (DM_DA*) dm->data; 95 const PetscInt dim = da->dim; 96 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 97 const PetscInt nVx = mx+1; 98 const PetscInt nVy = dim > 1 ? (my+1) : 1; 99 const PetscInt nVz = dim > 2 ? (mz+1) : 1; 100 const PetscInt nV = nVx*nVy*nVz; 101 102 PetscFunctionBegin; 103 if (numVerticesX) { 104 PetscValidIntPointer(numVerticesX,2); 105 *numVerticesX = nVx; 106 } 107 if (numVerticesY) { 108 PetscValidIntPointer(numVerticesY,3); 109 *numVerticesY = nVy; 110 } 111 if (numVerticesZ) { 112 PetscValidIntPointer(numVerticesZ,4); 113 *numVerticesZ = nVz; 114 } 115 if (numVertices) { 116 PetscValidIntPointer(numVertices,5); 117 *numVertices = nV; 118 } 119 PetscFunctionReturn(0); 120 } 121 122 #undef __FUNCT__ 123 #define __FUNCT__ "DMDAGetNumFaces" 124 PetscErrorCode DMDAGetNumFaces(DM dm, PetscInt *numXFacesX, PetscInt *numXFaces, PetscInt *numYFacesY, PetscInt *numYFaces, PetscInt *numZFacesZ, PetscInt *numZFaces) 125 { 126 DM_DA *da = (DM_DA*) dm->data; 127 const PetscInt dim = da->dim; 128 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 129 const PetscInt nxF = (dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1); 130 const PetscInt nXF = (mx+1)*nxF; 131 const PetscInt nyF = mx*(dim > 2 ? mz : 1); 132 const PetscInt nYF = dim > 1 ? (my+1)*nyF : 0; 133 const PetscInt nzF = mx*(dim > 1 ? my : 0); 134 const PetscInt nZF = dim > 2 ? (mz+1)*nzF : 0; 135 136 PetscFunctionBegin; 137 if (numXFacesX) { 138 PetscValidIntPointer(numXFacesX,2); 139 *numXFacesX = nxF; 140 } 141 if (numXFaces) { 142 PetscValidIntPointer(numXFaces,3); 143 *numXFaces = nXF; 144 } 145 if (numYFacesY) { 146 PetscValidIntPointer(numYFacesY,4); 147 *numYFacesY = nyF; 148 } 149 if (numYFaces) { 150 PetscValidIntPointer(numYFaces,5); 151 *numYFaces = nYF; 152 } 153 if (numZFacesZ) { 154 PetscValidIntPointer(numZFacesZ,6); 155 *numZFacesZ = nzF; 156 } 157 if (numZFaces) { 158 PetscValidIntPointer(numZFaces,7); 159 *numZFaces = nZF; 160 } 161 PetscFunctionReturn(0); 162 } 163 164 #undef __FUNCT__ 165 #define __FUNCT__ "DMDAGetHeightStratum" 166 PetscErrorCode DMDAGetHeightStratum(DM dm, PetscInt height, PetscInt *pStart, PetscInt *pEnd) 167 { 168 DM_DA *da = (DM_DA*) dm->data; 169 const PetscInt dim = da->dim; 170 PetscInt nC, nV, nXF, nYF, nZF; 171 PetscErrorCode ierr; 172 173 PetscFunctionBegin; 174 if (pStart) PetscValidIntPointer(pStart,3); 175 if (pEnd) PetscValidIntPointer(pEnd,4); 176 ierr = DMDAGetNumCells(dm, &nC);CHKERRQ(ierr); 177 ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr); 178 ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr); 179 if (height == 0) { 180 /* Cells */ 181 if (pStart) *pStart = 0; 182 if (pEnd) *pEnd = nC; 183 } else if (height == 1) { 184 /* Faces */ 185 if (pStart) *pStart = nC+nV; 186 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 187 } else if (height == dim) { 188 /* Vertices */ 189 if (pStart) *pStart = nC; 190 if (pEnd) *pEnd = nC+nV; 191 } else if (height < 0) { 192 /* All points */ 193 if (pStart) *pStart = 0; 194 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 195 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of height %d in the DA", height); 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 NULL for 1 209 . numVertexDof - The number of dofs per vertex for each field, or NULL 210 . numFaceDof - The number of dofs per face for each field and direction, or NULL 211 - numCellDof - The number of dofs per cell for each field, or 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(PetscObjectComm((PetscObject)dm), &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(PetscObjectComm((PetscObject)dm), 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) { 262 numFaceTotDof[0] += numFaceDof[f*dim+0]; 263 numFaceTotDof[1] += dim > 1 ? numFaceDof[f*dim+1] : 0; 264 numFaceTotDof[2] += dim > 2 ? numFaceDof[f*dim+2] : 0; 265 } 266 } 267 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &dm->defaultSection);CHKERRQ(ierr); 268 if (numFields > 1) { 269 ierr = PetscSectionSetNumFields(dm->defaultSection, numFields);CHKERRQ(ierr); 270 for (f = 0; f < numFields; ++f) { 271 const char *name; 272 273 ierr = DMDAGetFieldName(dm, f, &name);CHKERRQ(ierr); 274 ierr = PetscSectionSetFieldName(dm->defaultSection, f, name);CHKERRQ(ierr); 275 if (numComp) { 276 ierr = PetscSectionSetFieldComponents(dm->defaultSection, f, numComp[f]);CHKERRQ(ierr); 277 } 278 } 279 } else { 280 numFields = 0; 281 } 282 ierr = PetscSectionSetChart(dm->defaultSection, pStart, pEnd);CHKERRQ(ierr); 283 if (numVertexDof) { 284 for (v = vStart; v < vEnd; ++v) { 285 for (f = 0; f < numFields; ++f) { 286 ierr = PetscSectionSetFieldDof(dm->defaultSection, v, f, numVertexDof[f]);CHKERRQ(ierr); 287 } 288 ierr = PetscSectionSetDof(dm->defaultSection, v, numVertexTotDof);CHKERRQ(ierr); 289 } 290 } 291 if (numFaceDof) { 292 for (xf = xfStart; xf < xfEnd; ++xf) { 293 for (f = 0; f < numFields; ++f) { 294 ierr = PetscSectionSetFieldDof(dm->defaultSection, xf, f, numFaceDof[f*dim+0]);CHKERRQ(ierr); 295 } 296 ierr = PetscSectionSetDof(dm->defaultSection, xf, numFaceTotDof[0]);CHKERRQ(ierr); 297 } 298 for (yf = yfStart; yf < yfEnd; ++yf) { 299 for (f = 0; f < numFields; ++f) { 300 ierr = PetscSectionSetFieldDof(dm->defaultSection, yf, f, numFaceDof[f*dim+1]);CHKERRQ(ierr); 301 } 302 ierr = PetscSectionSetDof(dm->defaultSection, yf, numFaceTotDof[1]);CHKERRQ(ierr); 303 } 304 for (zf = zfStart; zf < zfEnd; ++zf) { 305 for (f = 0; f < numFields; ++f) { 306 ierr = PetscSectionSetFieldDof(dm->defaultSection, zf, f, numFaceDof[f*dim+2]);CHKERRQ(ierr); 307 } 308 ierr = PetscSectionSetDof(dm->defaultSection, zf, numFaceTotDof[2]);CHKERRQ(ierr); 309 } 310 } 311 if (numCellDof) { 312 for (c = cStart; c < cEnd; ++c) { 313 for (f = 0; f < numFields; ++f) { 314 ierr = PetscSectionSetFieldDof(dm->defaultSection, c, f, numCellDof[f]);CHKERRQ(ierr); 315 } 316 ierr = PetscSectionSetDof(dm->defaultSection, c, numCellTotDof);CHKERRQ(ierr); 317 } 318 } 319 ierr = PetscSectionSetUp(dm->defaultSection);CHKERRQ(ierr); 320 /* Create mesh point SF */ 321 ierr = DMDAGetNeighbors(dm, &neighbors);CHKERRQ(ierr); 322 for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) { 323 for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) { 324 for (xn = 0; xn < 3; ++xn) { 325 const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0; 326 const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn]; 327 328 if (neighbor >= 0 && neighbor < rank) { 329 nleaves += (!xp ? nVx : 1) * (!yp ? nVy : 1) * (!zp ? nVz : 1); /* vertices */ 330 if (xp && !yp && !zp) { 331 nleaves += nxF; /* x faces */ 332 } else if (yp && !zp && !xp) { 333 nleaves += nyF; /* y faces */ 334 } else if (zp && !xp && !yp) { 335 nleaves += nzF; /* z faces */ 336 } 337 } 338 } 339 } 340 } 341 ierr = PetscMalloc2(nleaves,PetscInt,&localPoints,nleaves,PetscSFNode,&remotePoints);CHKERRQ(ierr); 342 for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) { 343 for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) { 344 for (xn = 0; xn < 3; ++xn) { 345 const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0; 346 const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn]; 347 PetscInt xv, yv, zv; 348 349 if (neighbor >= 0 && neighbor < rank) { 350 if (xp < 0) { /* left */ 351 if (yp < 0) { /* bottom */ 352 if (zp < 0) { /* back */ 353 const PetscInt localVertex = ( 0*nVy + 0)*nVx + 0 + nC; 354 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 355 nleavesCheck += 1; /* left bottom back vertex */ 356 357 localPoints[nL] = localVertex; 358 remotePoints[nL].rank = neighbor; 359 remotePoints[nL].index = remoteVertex; 360 ++nL; 361 } else if (zp > 0) { /* front */ 362 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; 363 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 364 nleavesCheck += 1; /* left bottom front vertex */ 365 366 localPoints[nL] = localVertex; 367 remotePoints[nL].rank = neighbor; 368 remotePoints[nL].index = remoteVertex; 369 ++nL; 370 } else { 371 nleavesCheck += nVz; /* left bottom vertices */ 372 for (zv = 0; zv < nVz; ++zv, ++nL) { 373 const PetscInt localVertex = (zv*nVy + 0)*nVx + 0 + nC; 374 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 375 376 localPoints[nL] = localVertex; 377 remotePoints[nL].rank = neighbor; 378 remotePoints[nL].index = remoteVertex; 379 } 380 } 381 } else if (yp > 0) { /* top */ 382 if (zp < 0) { /* back */ 383 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; 384 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 385 nleavesCheck += 1; /* left top back vertex */ 386 387 localPoints[nL] = localVertex; 388 remotePoints[nL].rank = neighbor; 389 remotePoints[nL].index = remoteVertex; 390 ++nL; 391 } else if (zp > 0) { /* front */ 392 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; 393 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 394 nleavesCheck += 1; /* left top front vertex */ 395 396 localPoints[nL] = localVertex; 397 remotePoints[nL].rank = neighbor; 398 remotePoints[nL].index = remoteVertex; 399 ++nL; 400 } else { 401 nleavesCheck += nVz; /* left top vertices */ 402 for (zv = 0; zv < nVz; ++zv, ++nL) { 403 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; 404 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 405 406 localPoints[nL] = localVertex; 407 remotePoints[nL].rank = neighbor; 408 remotePoints[nL].index = remoteVertex; 409 } 410 } 411 } else { 412 if (zp < 0) { /* back */ 413 nleavesCheck += nVy; /* left back vertices */ 414 for (yv = 0; yv < nVy; ++yv, ++nL) { 415 const PetscInt localVertex = ( 0*nVy + yv)*nVx + 0 + nC; 416 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 417 418 localPoints[nL] = localVertex; 419 remotePoints[nL].rank = neighbor; 420 remotePoints[nL].index = remoteVertex; 421 } 422 } else if (zp > 0) { /* front */ 423 nleavesCheck += nVy; /* left front vertices */ 424 for (yv = 0; yv < nVy; ++yv, ++nL) { 425 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; 426 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 427 428 localPoints[nL] = localVertex; 429 remotePoints[nL].rank = neighbor; 430 remotePoints[nL].index = remoteVertex; 431 } 432 } else { 433 nleavesCheck += nVy*nVz; /* left vertices */ 434 for (zv = 0; zv < nVz; ++zv) { 435 for (yv = 0; yv < nVy; ++yv, ++nL) { 436 const PetscInt localVertex = (zv*nVy + yv)*nVx + 0 + nC; 437 const PetscInt remoteVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 438 439 localPoints[nL] = localVertex; 440 remotePoints[nL].rank = neighbor; 441 remotePoints[nL].index = remoteVertex; 442 } 443 } 444 nleavesCheck += nxF; /* left faces */ 445 for (xf = 0; xf < nxF; ++xf, ++nL) { 446 /* THIS IS WRONG */ 447 const PetscInt localFace = 0 + nC+nV; 448 const PetscInt remoteFace = 0 + nC+nV; 449 450 localPoints[nL] = localFace; 451 remotePoints[nL].rank = neighbor; 452 remotePoints[nL].index = remoteFace; 453 } 454 } 455 } 456 } else if (xp > 0) { /* right */ 457 if (yp < 0) { /* bottom */ 458 if (zp < 0) { /* back */ 459 const PetscInt localVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; 460 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 461 nleavesCheck += 1; /* right bottom back vertex */ 462 463 localPoints[nL] = localVertex; 464 remotePoints[nL].rank = neighbor; 465 remotePoints[nL].index = remoteVertex; 466 ++nL; 467 } else if (zp > 0) { /* front */ 468 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; 469 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 470 nleavesCheck += 1; /* right bottom front vertex */ 471 472 localPoints[nL] = localVertex; 473 remotePoints[nL].rank = neighbor; 474 remotePoints[nL].index = remoteVertex; 475 ++nL; 476 } else { 477 nleavesCheck += nVz; /* right bottom vertices */ 478 for (zv = 0; zv < nVz; ++zv, ++nL) { 479 const PetscInt localVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; 480 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 481 482 localPoints[nL] = localVertex; 483 remotePoints[nL].rank = neighbor; 484 remotePoints[nL].index = remoteVertex; 485 } 486 } 487 } else if (yp > 0) { /* top */ 488 if (zp < 0) { /* back */ 489 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; 490 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 491 nleavesCheck += 1; /* right top back vertex */ 492 493 localPoints[nL] = localVertex; 494 remotePoints[nL].rank = neighbor; 495 remotePoints[nL].index = remoteVertex; 496 ++nL; 497 } else if (zp > 0) { /* front */ 498 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; 499 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 500 nleavesCheck += 1; /* right top front vertex */ 501 502 localPoints[nL] = localVertex; 503 remotePoints[nL].rank = neighbor; 504 remotePoints[nL].index = remoteVertex; 505 ++nL; 506 } else { 507 nleavesCheck += nVz; /* right top vertices */ 508 for (zv = 0; zv < nVz; ++zv, ++nL) { 509 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; 510 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 511 512 localPoints[nL] = localVertex; 513 remotePoints[nL].rank = neighbor; 514 remotePoints[nL].index = remoteVertex; 515 } 516 } 517 } else { 518 if (zp < 0) { /* back */ 519 nleavesCheck += nVy; /* right back vertices */ 520 for (yv = 0; yv < nVy; ++yv, ++nL) { 521 const PetscInt localVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; 522 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 523 524 localPoints[nL] = localVertex; 525 remotePoints[nL].rank = neighbor; 526 remotePoints[nL].index = remoteVertex; 527 } 528 } else if (zp > 0) { /* front */ 529 nleavesCheck += nVy; /* right front vertices */ 530 for (yv = 0; yv < nVy; ++yv, ++nL) { 531 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; 532 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 533 534 localPoints[nL] = localVertex; 535 remotePoints[nL].rank = neighbor; 536 remotePoints[nL].index = remoteVertex; 537 } 538 } else { 539 nleavesCheck += nVy*nVz; /* right vertices */ 540 for (zv = 0; zv < nVz; ++zv) { 541 for (yv = 0; yv < nVy; ++yv, ++nL) { 542 const PetscInt localVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; 543 const PetscInt remoteVertex = (zv*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 544 545 localPoints[nL] = localVertex; 546 remotePoints[nL].rank = neighbor; 547 remotePoints[nL].index = remoteVertex; 548 } 549 } 550 nleavesCheck += nxF; /* right faces */ 551 for (xf = 0; xf < nxF; ++xf, ++nL) { 552 /* THIS IS WRONG */ 553 const PetscInt localFace = 0 + nC+nV; 554 const PetscInt remoteFace = 0 + nC+nV; 555 556 localPoints[nL] = localFace; 557 remotePoints[nL].rank = neighbor; 558 remotePoints[nL].index = remoteFace; 559 } 560 } 561 } 562 } else { 563 if (yp < 0) { /* bottom */ 564 if (zp < 0) { /* back */ 565 nleavesCheck += nVx; /* bottom back vertices */ 566 for (xv = 0; xv < nVx; ++xv, ++nL) { 567 const PetscInt localVertex = ( 0*nVy + 0)*nVx + xv + nC; 568 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 569 570 localPoints[nL] = localVertex; 571 remotePoints[nL].rank = neighbor; 572 remotePoints[nL].index = remoteVertex; 573 } 574 } else if (zp > 0) { /* front */ 575 nleavesCheck += nVx; /* bottom front vertices */ 576 for (xv = 0; xv < nVx; ++xv, ++nL) { 577 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; 578 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 579 580 localPoints[nL] = localVertex; 581 remotePoints[nL].rank = neighbor; 582 remotePoints[nL].index = remoteVertex; 583 } 584 } else { 585 nleavesCheck += nVx*nVz; /* bottom vertices */ 586 for (zv = 0; zv < nVz; ++zv) { 587 for (xv = 0; xv < nVx; ++xv, ++nL) { 588 const PetscInt localVertex = (zv*nVy + 0)*nVx + xv + nC; 589 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 590 591 localPoints[nL] = localVertex; 592 remotePoints[nL].rank = neighbor; 593 remotePoints[nL].index = remoteVertex; 594 } 595 } 596 nleavesCheck += nyF; /* bottom faces */ 597 for (yf = 0; yf < nyF; ++yf, ++nL) { 598 /* THIS IS WRONG */ 599 const PetscInt localFace = 0 + nC+nV; 600 const PetscInt remoteFace = 0 + nC+nV; 601 602 localPoints[nL] = localFace; 603 remotePoints[nL].rank = neighbor; 604 remotePoints[nL].index = remoteFace; 605 } 606 } 607 } else if (yp > 0) { /* top */ 608 if (zp < 0) { /* back */ 609 nleavesCheck += nVx; /* top back vertices */ 610 for (xv = 0; xv < nVx; ++xv, ++nL) { 611 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; 612 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 613 614 localPoints[nL] = localVertex; 615 remotePoints[nL].rank = neighbor; 616 remotePoints[nL].index = remoteVertex; 617 } 618 } else if (zp > 0) { /* front */ 619 nleavesCheck += nVx; /* top front vertices */ 620 for (xv = 0; xv < nVx; ++xv, ++nL) { 621 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; 622 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 623 624 localPoints[nL] = localVertex; 625 remotePoints[nL].rank = neighbor; 626 remotePoints[nL].index = remoteVertex; 627 } 628 } else { 629 nleavesCheck += nVx*nVz; /* top vertices */ 630 for (zv = 0; zv < nVz; ++zv) { 631 for (xv = 0; xv < nVx; ++xv, ++nL) { 632 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + xv + nC; 633 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 634 635 localPoints[nL] = localVertex; 636 remotePoints[nL].rank = neighbor; 637 remotePoints[nL].index = remoteVertex; 638 } 639 } 640 nleavesCheck += nyF; /* top faces */ 641 for (yf = 0; yf < nyF; ++yf, ++nL) { 642 /* THIS IS WRONG */ 643 const PetscInt localFace = 0 + nC+nV; 644 const PetscInt remoteFace = 0 + nC+nV; 645 646 localPoints[nL] = localFace; 647 remotePoints[nL].rank = neighbor; 648 remotePoints[nL].index = remoteFace; 649 } 650 } 651 } else { 652 if (zp < 0) { /* back */ 653 nleavesCheck += nVx*nVy; /* back vertices */ 654 for (yv = 0; yv < nVy; ++yv) { 655 for (xv = 0; xv < nVx; ++xv, ++nL) { 656 const PetscInt localVertex = ( 0*nVy + yv)*nVx + xv + nC; 657 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 658 659 localPoints[nL] = localVertex; 660 remotePoints[nL].rank = neighbor; 661 remotePoints[nL].index = remoteVertex; 662 } 663 } 664 nleavesCheck += nzF; /* back faces */ 665 for (zf = 0; zf < nzF; ++zf, ++nL) { 666 /* THIS IS WRONG */ 667 const PetscInt localFace = 0 + nC+nV; 668 const PetscInt remoteFace = 0 + nC+nV; 669 670 localPoints[nL] = localFace; 671 remotePoints[nL].rank = neighbor; 672 remotePoints[nL].index = remoteFace; 673 } 674 } else if (zp > 0) { /* front */ 675 nleavesCheck += nVx*nVy; /* front vertices */ 676 for (yv = 0; yv < nVy; ++yv) { 677 for (xv = 0; xv < nVx; ++xv, ++nL) { 678 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; 679 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 680 681 localPoints[nL] = localVertex; 682 remotePoints[nL].rank = neighbor; 683 remotePoints[nL].index = remoteVertex; 684 } 685 } 686 nleavesCheck += nzF; /* front faces */ 687 for (zf = 0; zf < nzF; ++zf, ++nL) { 688 /* THIS IS WRONG */ 689 const PetscInt localFace = 0 + nC+nV; 690 const PetscInt remoteFace = 0 + nC+nV; 691 692 localPoints[nL] = localFace; 693 remotePoints[nL].rank = neighbor; 694 remotePoints[nL].index = remoteFace; 695 } 696 } else { 697 /* Nothing is shared from the interior */ 698 } 699 } 700 } 701 } 702 } 703 } 704 } 705 /* TODO: Remove duplication in leaf determination */ 706 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); 707 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf);CHKERRQ(ierr); 708 ierr = PetscSFSetGraph(sf, pEnd, nleaves, localPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 709 /* Create global section */ 710 ierr = PetscSectionCreateGlobalSection(dm->defaultSection, sf, PETSC_FALSE, &dm->defaultGlobalSection);CHKERRQ(ierr); 711 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 712 /* Create default SF */ 713 ierr = DMCreateDefaultSF(dm, dm->defaultSection, dm->defaultGlobalSection);CHKERRQ(ierr); 714 PetscFunctionReturn(0); 715 } 716 717 /* ------------------------------------------------------------------- */ 718 719 #undef __FUNCT__ 720 #define __FUNCT__ "DMDAGetArray" 721 /*@C 722 DMDAGetArray - Gets a work array for a DMDA 723 724 Input Parameter: 725 + da - information about my local patch 726 - ghosted - do you want arrays for the ghosted or nonghosted patch 727 728 Output Parameters: 729 . vptr - array data structured 730 731 Note: The vector values are NOT initialized and may have garbage in them, so you may need 732 to zero them. 733 734 Level: advanced 735 736 .seealso: DMDARestoreArray() 737 738 @*/ 739 PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr) 740 { 741 PetscErrorCode ierr; 742 PetscInt j,i,xs,ys,xm,ym,zs,zm; 743 char *iarray_start; 744 void **iptr = (void**)vptr; 745 DM_DA *dd = (DM_DA*)da->data; 746 747 PetscFunctionBegin; 748 PetscValidHeaderSpecific(da,DM_CLASSID,1); 749 if (ghosted) { 750 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 751 if (dd->arrayghostedin[i]) { 752 *iptr = dd->arrayghostedin[i]; 753 iarray_start = (char*)dd->startghostedin[i]; 754 dd->arrayghostedin[i] = NULL; 755 dd->startghostedin[i] = NULL; 756 757 goto done; 758 } 759 } 760 xs = dd->Xs; 761 ys = dd->Ys; 762 zs = dd->Zs; 763 xm = dd->Xe-dd->Xs; 764 ym = dd->Ye-dd->Ys; 765 zm = dd->Ze-dd->Zs; 766 } else { 767 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 768 if (dd->arrayin[i]) { 769 *iptr = dd->arrayin[i]; 770 iarray_start = (char*)dd->startin[i]; 771 dd->arrayin[i] = NULL; 772 dd->startin[i] = NULL; 773 774 goto done; 775 } 776 } 777 xs = dd->xs; 778 ys = dd->ys; 779 zs = dd->zs; 780 xm = dd->xe-dd->xs; 781 ym = dd->ye-dd->ys; 782 zm = dd->ze-dd->zs; 783 } 784 785 switch (dd->dim) { 786 case 1: { 787 void *ptr; 788 789 ierr = PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 790 791 ptr = (void*)(iarray_start - xs*sizeof(PetscScalar)); 792 *iptr = (void*)ptr; 793 break; 794 } 795 case 2: { 796 void **ptr; 797 798 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 799 800 ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*)); 801 for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs); 802 *iptr = (void*)ptr; 803 break; 804 } 805 case 3: { 806 void ***ptr,**bptr; 807 808 ierr = PetscMalloc((zm+1)*sizeof(void**)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 809 810 ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*)); 811 bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**)); 812 for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys); 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 *iptr = (void*)ptr; 819 break; 820 } 821 default: 822 SETERRQ1(PetscObjectComm((PetscObject)da),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] = NULL; 875 dd->startghostedout[i] = 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] = NULL; 891 dd->startout[i] = 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