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