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 <petscbt.h> 8 #include <petscsf.h> 9 10 /* 11 This allows the DMDA vectors to properly tell MATLAB their dimensions 12 */ 13 #if defined(PETSC_HAVE_MATLAB_ENGINE) 14 #include <engine.h> /* MATLAB include file */ 15 #include <mex.h> /* MATLAB include file */ 16 #undef __FUNCT__ 17 #define __FUNCT__ "VecMatlabEnginePut_DA2d" 18 static PetscErrorCode VecMatlabEnginePut_DA2d(PetscObject obj,void *mengine) 19 { 20 PetscErrorCode ierr; 21 PetscInt n,m; 22 Vec vec = (Vec)obj; 23 PetscScalar *array; 24 mxArray *mat; 25 DM da; 26 27 PetscFunctionBegin; 28 ierr = VecGetDM(vec, &da);CHKERRQ(ierr); 29 if (!da) SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_ARG_WRONGSTATE,"Vector not associated with a DMDA"); 30 ierr = DMDAGetGhostCorners(da,0,0,0,&m,&n,0);CHKERRQ(ierr); 31 32 ierr = VecGetArray(vec,&array);CHKERRQ(ierr); 33 #if !defined(PETSC_USE_COMPLEX) 34 mat = mxCreateDoubleMatrix(m,n,mxREAL); 35 #else 36 mat = mxCreateDoubleMatrix(m,n,mxCOMPLEX); 37 #endif 38 ierr = PetscMemcpy(mxGetPr(mat),array,n*m*sizeof(PetscScalar));CHKERRQ(ierr); 39 ierr = PetscObjectName(obj);CHKERRQ(ierr); 40 engPutVariable((Engine*)mengine,obj->name,mat); 41 42 ierr = VecRestoreArray(vec,&array);CHKERRQ(ierr); 43 PetscFunctionReturn(0); 44 } 45 #endif 46 47 48 #undef __FUNCT__ 49 #define __FUNCT__ "DMCreateLocalVector_DA" 50 PetscErrorCode DMCreateLocalVector_DA(DM da,Vec *g) 51 { 52 PetscErrorCode ierr; 53 DM_DA *dd = (DM_DA*)da->data; 54 55 PetscFunctionBegin; 56 PetscValidHeaderSpecific(da,DM_CLASSID,1); 57 PetscValidPointer(g,2); 58 if (da->defaultSection) { 59 ierr = DMCreateLocalVector_Section_Private(da,g);CHKERRQ(ierr); 60 } else { 61 ierr = VecCreate(PETSC_COMM_SELF,g);CHKERRQ(ierr); 62 ierr = VecSetSizes(*g,dd->nlocal,PETSC_DETERMINE);CHKERRQ(ierr); 63 ierr = VecSetBlockSize(*g,dd->w);CHKERRQ(ierr); 64 ierr = VecSetType(*g,da->vectype);CHKERRQ(ierr); 65 ierr = VecSetDM(*g, da);CHKERRQ(ierr); 66 #if defined(PETSC_HAVE_MATLAB_ENGINE) 67 if (dd->w == 1 && dd->dim == 2) { 68 ierr = PetscObjectComposeFunction((PetscObject)*g,"PetscMatlabEnginePut_C",VecMatlabEnginePut_DA2d);CHKERRQ(ierr); 69 } 70 #endif 71 } 72 PetscFunctionReturn(0); 73 } 74 75 #undef __FUNCT__ 76 #define __FUNCT__ "DMDAGetNumCells" 77 PetscErrorCode DMDAGetNumCells(DM dm, PetscInt *numCellsX, PetscInt *numCellsY, PetscInt *numCellsZ, PetscInt *numCells) 78 { 79 DM_DA *da = (DM_DA*) dm->data; 80 const PetscInt dim = da->dim; 81 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 82 const PetscInt nC = (mx)*(dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1); 83 84 PetscFunctionBegin; 85 if (numCellsX) { 86 PetscValidIntPointer(numCellsX,2); 87 *numCellsX = mx; 88 } 89 if (numCellsY) { 90 PetscValidIntPointer(numCellsX,3); 91 *numCellsY = my; 92 } 93 if (numCellsZ) { 94 PetscValidIntPointer(numCellsX,4); 95 *numCellsZ = mz; 96 } 97 if (numCells) { 98 PetscValidIntPointer(numCells,5); 99 *numCells = nC; 100 } 101 PetscFunctionReturn(0); 102 } 103 104 #undef __FUNCT__ 105 #define __FUNCT__ "DMDAGetNumVertices" 106 PetscErrorCode DMDAGetNumVertices(DM dm, PetscInt *numVerticesX, PetscInt *numVerticesY, PetscInt *numVerticesZ, PetscInt *numVertices) 107 { 108 DM_DA *da = (DM_DA*) dm->data; 109 const PetscInt dim = da->dim; 110 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 111 const PetscInt nVx = mx+1; 112 const PetscInt nVy = dim > 1 ? (my+1) : 1; 113 const PetscInt nVz = dim > 2 ? (mz+1) : 1; 114 const PetscInt nV = nVx*nVy*nVz; 115 116 PetscFunctionBegin; 117 if (numVerticesX) { 118 PetscValidIntPointer(numVerticesX,2); 119 *numVerticesX = nVx; 120 } 121 if (numVerticesY) { 122 PetscValidIntPointer(numVerticesY,3); 123 *numVerticesY = nVy; 124 } 125 if (numVerticesZ) { 126 PetscValidIntPointer(numVerticesZ,4); 127 *numVerticesZ = nVz; 128 } 129 if (numVertices) { 130 PetscValidIntPointer(numVertices,5); 131 *numVertices = nV; 132 } 133 PetscFunctionReturn(0); 134 } 135 136 #undef __FUNCT__ 137 #define __FUNCT__ "DMDAGetNumFaces" 138 PetscErrorCode DMDAGetNumFaces(DM dm, PetscInt *numXFacesX, PetscInt *numXFaces, PetscInt *numYFacesY, PetscInt *numYFaces, PetscInt *numZFacesZ, PetscInt *numZFaces) 139 { 140 DM_DA *da = (DM_DA*) dm->data; 141 const PetscInt dim = da->dim; 142 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 143 const PetscInt nxF = (dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1); 144 const PetscInt nXF = (mx+1)*nxF; 145 const PetscInt nyF = mx*(dim > 2 ? mz : 1); 146 const PetscInt nYF = dim > 1 ? (my+1)*nyF : 0; 147 const PetscInt nzF = mx*(dim > 1 ? my : 0); 148 const PetscInt nZF = dim > 2 ? (mz+1)*nzF : 0; 149 150 PetscFunctionBegin; 151 if (numXFacesX) { 152 PetscValidIntPointer(numXFacesX,2); 153 *numXFacesX = nxF; 154 } 155 if (numXFaces) { 156 PetscValidIntPointer(numXFaces,3); 157 *numXFaces = nXF; 158 } 159 if (numYFacesY) { 160 PetscValidIntPointer(numYFacesY,4); 161 *numYFacesY = nyF; 162 } 163 if (numYFaces) { 164 PetscValidIntPointer(numYFaces,5); 165 *numYFaces = nYF; 166 } 167 if (numZFacesZ) { 168 PetscValidIntPointer(numZFacesZ,6); 169 *numZFacesZ = nzF; 170 } 171 if (numZFaces) { 172 PetscValidIntPointer(numZFaces,7); 173 *numZFaces = nZF; 174 } 175 PetscFunctionReturn(0); 176 } 177 178 #undef __FUNCT__ 179 #define __FUNCT__ "DMDAGetHeightStratum" 180 PetscErrorCode DMDAGetHeightStratum(DM dm, PetscInt height, PetscInt *pStart, PetscInt *pEnd) 181 { 182 DM_DA *da = (DM_DA*) dm->data; 183 const PetscInt dim = da->dim; 184 PetscInt nC, nV, nXF, nYF, nZF; 185 PetscErrorCode ierr; 186 187 PetscFunctionBegin; 188 if (pStart) PetscValidIntPointer(pStart,3); 189 if (pEnd) PetscValidIntPointer(pEnd,4); 190 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); 191 ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr); 192 ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr); 193 if (height == 0) { 194 /* Cells */ 195 if (pStart) *pStart = 0; 196 if (pEnd) *pEnd = nC; 197 } else if (height == 1) { 198 /* Faces */ 199 if (pStart) *pStart = nC+nV; 200 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 201 } else if (height == dim) { 202 /* Vertices */ 203 if (pStart) *pStart = nC; 204 if (pEnd) *pEnd = nC+nV; 205 } else if (height < 0) { 206 /* All points */ 207 if (pStart) *pStart = 0; 208 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 209 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of height %d in the DA", height); 210 PetscFunctionReturn(0); 211 } 212 213 #undef __FUNCT__ 214 #define __FUNCT__ "DMDAGetDepthStratum" 215 PetscErrorCode DMDAGetDepthStratum(DM dm, PetscInt depth, PetscInt *pStart, PetscInt *pEnd) 216 { 217 DM_DA *da = (DM_DA*) dm->data; 218 const PetscInt dim = da->dim; 219 PetscInt nC, nV, nXF, nYF, nZF; 220 PetscErrorCode ierr; 221 222 PetscFunctionBegin; 223 if (pStart) PetscValidIntPointer(pStart,3); 224 if (pEnd) PetscValidIntPointer(pEnd,4); 225 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); 226 ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr); 227 ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr); 228 if (depth == dim) { 229 /* Cells */ 230 if (pStart) *pStart = 0; 231 if (pEnd) *pEnd = nC; 232 } else if (depth == dim-1) { 233 /* Faces */ 234 if (pStart) *pStart = nC+nV; 235 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 236 } else if (depth == 0) { 237 /* Vertices */ 238 if (pStart) *pStart = nC; 239 if (pEnd) *pEnd = nC+nV; 240 } else if (depth < 0) { 241 /* All points */ 242 if (pStart) *pStart = 0; 243 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 244 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of depth %d in the DA", depth); 245 PetscFunctionReturn(0); 246 } 247 248 #undef __FUNCT__ 249 #define __FUNCT__ "DMDAGetConeSize" 250 PetscErrorCode DMDAGetConeSize(DM dm, PetscInt p, PetscInt *coneSize) 251 { 252 DM_DA *da = (DM_DA*) dm->data; 253 const PetscInt dim = da->dim; 254 PetscInt nC, nV, nXF, nYF, nZF; 255 PetscErrorCode ierr; 256 257 PetscFunctionBegin; 258 *coneSize = 0; 259 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); 260 ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr); 261 ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr); 262 switch (dim) { 263 case 2: 264 if (p >= 0) { 265 if (p < nC) { 266 *coneSize = 4; 267 } else if (p < nC+nV) { 268 *coneSize = 0; 269 } else if (p < nC+nV+nXF+nYF+nZF) { 270 *coneSize = 2; 271 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF); 272 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p); 273 break; 274 case 3: 275 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D"); 276 break; 277 } 278 PetscFunctionReturn(0); 279 } 280 281 #undef __FUNCT__ 282 #define __FUNCT__ "DMDAGetCone" 283 PetscErrorCode DMDAGetCone(DM dm, PetscInt p, PetscInt *cone[]) 284 { 285 DM_DA *da = (DM_DA*) dm->data; 286 const PetscInt dim = da->dim; 287 PetscInt nCx, nCy, nCz, nC, nVx, nVy, nVz, nV, nxF, nyF, nzF, nXF, nYF, nZF; 288 PetscErrorCode ierr; 289 290 PetscFunctionBegin; 291 if (!cone) {ierr = DMGetWorkArray(dm, 6, PETSC_INT, cone);CHKERRQ(ierr);} 292 ierr = DMDAGetNumCells(dm, &nCx, &nCy, &nCz, &nC);CHKERRQ(ierr); 293 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr); 294 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 295 switch (dim) { 296 case 2: 297 if (p >= 0) { 298 if (p < nC) { 299 const PetscInt cy = p / nCx; 300 const PetscInt cx = p % nCx; 301 302 (*cone)[0] = cy*nxF + cx + nC+nV; 303 (*cone)[1] = cx*nyF + cy + nyF + nC+nV+nXF; 304 (*cone)[2] = cy*nxF + cx + nxF + nC+nV; 305 (*cone)[3] = cx*nyF + cy + nC+nV+nXF; 306 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do cell cones"); 307 } else if (p < nC+nV) { 308 } else if (p < nC+nV+nXF) { 309 const PetscInt fy = (p - nC+nV) / nxF; 310 const PetscInt fx = (p - nC+nV) % nxF; 311 312 (*cone)[0] = fy*nVx + fx + nC; 313 (*cone)[1] = fy*nVx + fx + 1 + nC; 314 } else if (p < nC+nV+nXF+nYF) { 315 const PetscInt fx = (p - nC+nV+nXF) / nyF; 316 const PetscInt fy = (p - nC+nV+nXF) % nyF; 317 318 (*cone)[0] = fy*nVx + fx + nC; 319 (*cone)[1] = fy*nVx + fx + nVx + nC; 320 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF); 321 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p); 322 break; 323 case 3: 324 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D"); 325 break; 326 } 327 PetscFunctionReturn(0); 328 } 329 330 #undef __FUNCT__ 331 #define __FUNCT__ "DMDARestoreCone" 332 PetscErrorCode DMDARestoreCone(DM dm, PetscInt p, PetscInt *cone[]) 333 { 334 PetscErrorCode ierr; 335 336 PetscFunctionBegin; 337 ierr = DMGetWorkArray(dm, 6, PETSC_INT, cone);CHKERRQ(ierr); 338 PetscFunctionReturn(0); 339 } 340 341 #undef __FUNCT__ 342 #define __FUNCT__ "DMDACreateSection" 343 /*@C 344 DMDACreateSection - Create a PetscSection inside the DMDA that describes data layout. This allows multiple fields with 345 different numbers of dofs on vertices, cells, and faces in each direction. 346 347 Input Parameters: 348 + dm- The DMDA 349 . numFields - The number of fields 350 . numComp - The number of components in each field 351 . numDof - The number of dofs per dimension for each field 352 . numFaceDof - The number of dofs per face for each field and direction, or NULL 353 354 Level: developer 355 356 Note: 357 The default DMDA numbering is as follows: 358 359 - Cells: [0, nC) 360 - Vertices: [nC, nC+nV) 361 - X-Faces: [nC+nV, nC+nV+nXF) normal is +- x-dir 362 - Y-Faces: [nC+nV+nXF, nC+nV+nXF+nYF) normal is +- y-dir 363 - Z-Faces: [nC+nV+nXF+nYF, nC+nV+nXF+nYF+nZF) normal is +- z-dir 364 365 We interpret the default DMDA partition as a cell partition, and the data assignment as a cell assignment. 366 @*/ 367 PetscErrorCode DMDACreateSection(DM dm, PetscInt numComp[], PetscInt numDof[], PetscInt numFaceDof[], PetscSection *s) 368 { 369 DM_DA *da = (DM_DA*) dm->data; 370 PetscSection section; 371 const PetscInt dim = da->dim; 372 PetscInt numFields, numVertexTotDof = 0, numCellTotDof = 0, numFaceTotDof[3] = {0, 0, 0}; 373 PetscBT isLeaf; 374 PetscSF sf; 375 PetscMPIInt rank; 376 const PetscMPIInt *neighbors; 377 PetscInt *localPoints; 378 PetscSFNode *remotePoints; 379 PetscInt nleaves = 0, nleavesCheck = 0, nL = 0; 380 PetscInt nC, nVx, nVy, nVz, nV, nxF, nXF, nyF, nYF, nzF, nZF; 381 PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart, zfEnd; 382 PetscInt f, v, c, xf, yf, zf, xn, yn, zn; 383 PetscErrorCode ierr; 384 385 PetscFunctionBegin; 386 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 387 PetscValidPointer(s, 4); 388 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 389 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); 390 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr); 391 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 392 ierr = DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);CHKERRQ(ierr); 393 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 394 ierr = DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 395 ierr = DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);CHKERRQ(ierr); 396 xfStart = vEnd; xfEnd = xfStart+nXF; 397 yfStart = xfEnd; yfEnd = yfStart+nYF; 398 zfStart = yfEnd; zfEnd = zfStart+nZF; 399 if (zfEnd != fEnd) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid face end %d, should be %d", zfEnd, fEnd); 400 /* Create local section */ 401 ierr = DMDAGetInfo(dm, 0,0,0,0,0,0,0, &numFields, 0,0,0,0,0);CHKERRQ(ierr); 402 for (f = 0; f < numFields; ++f) { 403 numVertexTotDof += numDof[f*(dim+1)+0]; 404 numCellTotDof += numDof[f*(dim+1)+dim]; 405 numFaceTotDof[0] += dim > 0 ? (numFaceDof ? numFaceDof[f*dim+0] : numDof[f*(dim+1)+dim-1]) : 0; 406 numFaceTotDof[1] += dim > 1 ? (numFaceDof ? numFaceDof[f*dim+1] : numDof[f*(dim+1)+dim-1]) : 0; 407 numFaceTotDof[2] += dim > 2 ? (numFaceDof ? numFaceDof[f*dim+2] : numDof[f*(dim+1)+dim-1]) : 0; 408 } 409 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion);CHKERRQ(ierr); 410 if (numFields > 0) { 411 ierr = PetscSectionSetNumFields(section, numFields);CHKERRQ(ierr); 412 for (f = 0; f < numFields; ++f) { 413 const char *name; 414 415 ierr = DMDAGetFieldName(dm, f, &name);CHKERRQ(ierr); 416 ierr = PetscSectionSetFieldName(section, f, name ? name : "Field");CHKERRQ(ierr); 417 if (numComp) { 418 ierr = PetscSectionSetFieldComponents(section, f, numComp[f]);CHKERRQ(ierr); 419 } 420 } 421 } 422 ierr = PetscSectionSetChart(section, pStart, pEnd);CHKERRQ(ierr); 423 for (v = vStart; v < vEnd; ++v) { 424 for (f = 0; f < numFields; ++f) { 425 ierr = PetscSectionSetFieldDof(section, v, f, numDof[f*(dim+1)+0]);CHKERRQ(ierr); 426 } 427 ierr = PetscSectionSetDof(section, v, numVertexTotDof);CHKERRQ(ierr); 428 } 429 for (xf = xfStart; xf < xfEnd; ++xf) { 430 for (f = 0; f < numFields; ++f) { 431 ierr = PetscSectionSetFieldDof(section, xf, f, numFaceDof ? numFaceDof[f*dim+0] : numDof[f*(dim+1)+dim-1]);CHKERRQ(ierr); 432 } 433 ierr = PetscSectionSetDof(section, xf, numFaceTotDof[0]);CHKERRQ(ierr); 434 } 435 for (yf = yfStart; yf < yfEnd; ++yf) { 436 for (f = 0; f < numFields; ++f) { 437 ierr = PetscSectionSetFieldDof(section, yf, f, numFaceDof ? numFaceDof[f*dim+1] : numDof[f*(dim+1)+dim-1]);CHKERRQ(ierr); 438 } 439 ierr = PetscSectionSetDof(section, yf, numFaceTotDof[1]);CHKERRQ(ierr); 440 } 441 for (zf = zfStart; zf < zfEnd; ++zf) { 442 for (f = 0; f < numFields; ++f) { 443 ierr = PetscSectionSetFieldDof(section, zf, f, numFaceDof ? numFaceDof[f*dim+2] : numDof[f*(dim+1)+dim-1]);CHKERRQ(ierr); 444 } 445 ierr = PetscSectionSetDof(section, zf, numFaceTotDof[2]);CHKERRQ(ierr); 446 } 447 for (c = cStart; c < cEnd; ++c) { 448 for (f = 0; f < numFields; ++f) { 449 ierr = PetscSectionSetFieldDof(section, c, f, numDof[f*(dim+1)+dim]);CHKERRQ(ierr); 450 } 451 ierr = PetscSectionSetDof(section, c, numCellTotDof);CHKERRQ(ierr); 452 } 453 ierr = PetscSectionSetUp(section);CHKERRQ(ierr); 454 /* Create mesh point SF */ 455 ierr = PetscBTCreate(pEnd-pStart, &isLeaf);CHKERRQ(ierr); 456 ierr = DMDAGetNeighbors(dm, &neighbors);CHKERRQ(ierr); 457 for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) { 458 for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) { 459 for (xn = 0; xn < 3; ++xn) { 460 const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0; 461 const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn]; 462 PetscInt xv, yv, zv; 463 464 if (neighbor >= 0 && neighbor < rank) { 465 if (xp < 0) { /* left */ 466 if (yp < 0) { /* bottom */ 467 if (zp < 0) { /* back */ 468 const PetscInt localVertex = ( 0*nVy + 0)*nVx + 0 + nC; /* left bottom back vertex */ 469 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 470 } else if (zp > 0) { /* front */ 471 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; /* left bottom front vertex */ 472 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 473 } else { 474 for (zv = 0; zv < nVz; ++zv) { 475 const PetscInt localVertex = (zv*nVy + 0)*nVx + 0 + nC; /* left bottom vertices */ 476 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 477 } 478 } 479 } else if (yp > 0) { /* top */ 480 if (zp < 0) { /* back */ 481 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; /* left top back vertex */ 482 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 483 } else if (zp > 0) { /* front */ 484 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; /* left top front vertex */ 485 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 486 } else { 487 for (zv = 0; zv < nVz; ++zv) { 488 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; /* left top vertices */ 489 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 490 } 491 } 492 } else { 493 if (zp < 0) { /* back */ 494 for (yv = 0; yv < nVy; ++yv) { 495 const PetscInt localVertex = ( 0*nVy + yv)*nVx + 0 + nC; /* left back vertices */ 496 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 497 } 498 } else if (zp > 0) { /* front */ 499 for (yv = 0; yv < nVy; ++yv) { 500 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; /* left front vertices */ 501 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 502 } 503 } else { 504 for (zv = 0; zv < nVz; ++zv) { 505 for (yv = 0; yv < nVy; ++yv) { 506 const PetscInt localVertex = (zv*nVy + yv)*nVx + 0 + nC; /* left vertices */ 507 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 508 } 509 } 510 #if 0 511 for (xf = 0; xf < nxF; ++xf) { 512 /* THIS IS WRONG */ 513 const PetscInt localFace = 0 + nC+nV; /* left faces */ 514 if (!PetscBTLookupSet(isLeaf, localFace)) ++nleaves; 515 } 516 #endif 517 } 518 } 519 } else if (xp > 0) { /* right */ 520 if (yp < 0) { /* bottom */ 521 if (zp < 0) { /* back */ 522 const PetscInt localVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; /* right bottom back vertex */ 523 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 524 } else if (zp > 0) { /* front */ 525 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; /* right bottom front vertex */ 526 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 527 } else { 528 for (zv = 0; zv < nVz; ++zv) { 529 const PetscInt localVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; /* right bottom vertices */ 530 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 531 } 532 } 533 } else if (yp > 0) { /* top */ 534 if (zp < 0) { /* back */ 535 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top back vertex */ 536 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 537 } else if (zp > 0) { /* front */ 538 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top front vertex */ 539 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 540 } else { 541 for (zv = 0; zv < nVz; ++zv) { 542 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top vertices */ 543 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 544 } 545 } 546 } else { 547 if (zp < 0) { /* back */ 548 for (yv = 0; yv < nVy; ++yv) { 549 const PetscInt localVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; /* right back vertices */ 550 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 551 } 552 } else if (zp > 0) { /* front */ 553 for (yv = 0; yv < nVy; ++yv) { 554 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* right front vertices */ 555 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 556 } 557 } else { 558 for (zv = 0; zv < nVz; ++zv) { 559 for (yv = 0; yv < nVy; ++yv) { 560 const PetscInt localVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* right vertices */ 561 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 562 } 563 } 564 #if 0 565 for (xf = 0; xf < nxF; ++xf) { 566 /* THIS IS WRONG */ 567 const PetscInt localFace = 0 + nC+nV; /* right faces */ 568 if (!PetscBTLookupSet(isLeaf, localFace)) ++nleaves; 569 } 570 #endif 571 } 572 } 573 } else { 574 if (yp < 0) { /* bottom */ 575 if (zp < 0) { /* back */ 576 for (xv = 0; xv < nVx; ++xv) { 577 const PetscInt localVertex = ( 0*nVy + 0)*nVx + xv + nC; /* bottom back vertices */ 578 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 579 } 580 } else if (zp > 0) { /* front */ 581 for (xv = 0; xv < nVx; ++xv) { 582 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; /* bottom front vertices */ 583 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 584 } 585 } else { 586 for (zv = 0; zv < nVz; ++zv) { 587 for (xv = 0; xv < nVx; ++xv) { 588 const PetscInt localVertex = (zv*nVy + 0)*nVx + xv + nC; /* bottom vertices */ 589 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 590 } 591 } 592 #if 0 593 for (yf = 0; yf < nyF; ++yf) { 594 /* THIS IS WRONG */ 595 const PetscInt localFace = 0 + nC+nV; /* bottom faces */ 596 if (!PetscBTLookupSet(isLeaf, localVFace)) ++nleaves; 597 } 598 #endif 599 } 600 } else if (yp > 0) { /* top */ 601 if (zp < 0) { /* back */ 602 for (xv = 0; xv < nVx; ++xv) { 603 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; /* top back vertices */ 604 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 605 } 606 } else if (zp > 0) { /* front */ 607 for (xv = 0; xv < nVx; ++xv) { 608 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* top front vertices */ 609 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 610 } 611 } else { 612 for (zv = 0; zv < nVz; ++zv) { 613 for (xv = 0; xv < nVx; ++xv) { 614 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* top vertices */ 615 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 616 } 617 } 618 #if 0 619 for (yf = 0; yf < nyF; ++yf) { 620 /* THIS IS WRONG */ 621 const PetscInt localFace = 0 + nC+nV; /* top faces */ 622 if (!PetscBTLookupSet(isLeaf, localVFace)) ++nleaves; 623 } 624 #endif 625 } 626 } else { 627 if (zp < 0) { /* back */ 628 for (yv = 0; yv < nVy; ++yv) { 629 for (xv = 0; xv < nVx; ++xv) { 630 const PetscInt localVertex = ( 0*nVy + yv)*nVx + xv + nC; /* back vertices */ 631 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 632 } 633 } 634 #if 0 635 for (zf = 0; zf < nzF; ++zf) { 636 /* THIS IS WRONG */ 637 const PetscInt localFace = 0 + nC+nV; /* back faces */ 638 if (!PetscBTLookupSet(isLeaf, localFace)) ++nleaves; 639 } 640 #endif 641 } else if (zp > 0) { /* front */ 642 for (yv = 0; yv < nVy; ++yv) { 643 for (xv = 0; xv < nVx; ++xv) { 644 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* front vertices */ 645 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves; 646 } 647 } 648 #if 0 649 for (zf = 0; zf < nzF; ++zf) { 650 /* THIS IS WRONG */ 651 const PetscInt localFace = 0 + nC+nV; /* front faces */ 652 if (!PetscBTLookupSet(isLeaf, localFace)) ++nleaves; 653 } 654 #endif 655 } else { 656 /* Nothing is shared from the interior */ 657 } 658 } 659 } 660 } 661 } 662 } 663 } 664 ierr = PetscBTMemzero(pEnd-pStart, isLeaf);CHKERRQ(ierr); 665 ierr = PetscMalloc2(nleaves,PetscInt,&localPoints,nleaves,PetscSFNode,&remotePoints);CHKERRQ(ierr); 666 for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) { 667 for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) { 668 for (xn = 0; xn < 3; ++xn) { 669 const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0; 670 const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn]; 671 PetscInt xv, yv, zv; 672 673 if (neighbor >= 0 && neighbor < rank) { 674 if (xp < 0) { /* left */ 675 if (yp < 0) { /* bottom */ 676 if (zp < 0) { /* back */ 677 const PetscInt localVertex = ( 0*nVy + 0)*nVx + 0 + nC; /* left bottom back vertex */ 678 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 679 680 if (!PetscBTLookupSet(isLeaf, localVertex)) { 681 localPoints[nL] = localVertex; 682 remotePoints[nL].rank = neighbor; 683 remotePoints[nL].index = remoteVertex; 684 ++nL; 685 } 686 } else if (zp > 0) { /* front */ 687 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; /* left bottom front vertex */ 688 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 689 690 if (!PetscBTLookupSet(isLeaf, localVertex)) { 691 localPoints[nL] = localVertex; 692 remotePoints[nL].rank = neighbor; 693 remotePoints[nL].index = remoteVertex; 694 ++nL; 695 } 696 } else { 697 for (zv = 0; zv < nVz; ++zv) { 698 const PetscInt localVertex = (zv*nVy + 0)*nVx + 0 + nC; /* left bottom vertices */ 699 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 700 701 if (!PetscBTLookupSet(isLeaf, localVertex)) { 702 localPoints[nL] = localVertex; 703 remotePoints[nL].rank = neighbor; 704 remotePoints[nL].index = remoteVertex; 705 ++nL; 706 } 707 } 708 } 709 } else if (yp > 0) { /* top */ 710 if (zp < 0) { /* back */ 711 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; /* left top back vertex */ 712 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 713 714 if (!PetscBTLookupSet(isLeaf, localVertex)) { 715 localPoints[nL] = localVertex; 716 remotePoints[nL].rank = neighbor; 717 remotePoints[nL].index = remoteVertex; 718 ++nL; 719 } 720 } else if (zp > 0) { /* front */ 721 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; /* left top front vertex */ 722 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 723 724 if (!PetscBTLookupSet(isLeaf, localVertex)) { 725 localPoints[nL] = localVertex; 726 remotePoints[nL].rank = neighbor; 727 remotePoints[nL].index = remoteVertex; 728 ++nL; 729 } 730 } else { 731 for (zv = 0; zv < nVz; ++zv) { 732 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; /* left top vertices */ 733 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 734 735 if (!PetscBTLookupSet(isLeaf, localVertex)) { 736 localPoints[nL] = localVertex; 737 remotePoints[nL].rank = neighbor; 738 remotePoints[nL].index = remoteVertex; 739 ++nL; 740 } 741 } 742 } 743 } else { 744 if (zp < 0) { /* back */ 745 for (yv = 0; yv < nVy; ++yv) { 746 const PetscInt localVertex = ( 0*nVy + yv)*nVx + 0 + nC; /* left back vertices */ 747 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 748 749 if (!PetscBTLookupSet(isLeaf, localVertex)) { 750 localPoints[nL] = localVertex; 751 remotePoints[nL].rank = neighbor; 752 remotePoints[nL].index = remoteVertex; 753 ++nL; 754 } 755 } 756 } else if (zp > 0) { /* front */ 757 for (yv = 0; yv < nVy; ++yv) { 758 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; /* left front vertices */ 759 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 760 761 if (!PetscBTLookupSet(isLeaf, localVertex)) { 762 localPoints[nL] = localVertex; 763 remotePoints[nL].rank = neighbor; 764 remotePoints[nL].index = remoteVertex; 765 ++nL; 766 } 767 } 768 } else { 769 for (zv = 0; zv < nVz; ++zv) { 770 for (yv = 0; yv < nVy; ++yv) { 771 const PetscInt localVertex = (zv*nVy + yv)*nVx + 0 + nC; /* left vertices */ 772 const PetscInt remoteVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 773 774 if (!PetscBTLookupSet(isLeaf, localVertex)) { 775 localPoints[nL] = localVertex; 776 remotePoints[nL].rank = neighbor; 777 remotePoints[nL].index = remoteVertex; 778 ++nL; 779 } 780 } 781 } 782 #if 0 783 for (xf = 0; xf < nxF; ++xf) { 784 /* THIS IS WRONG */ 785 const PetscInt localFace = 0 + nC+nV; /* left faces */ 786 const PetscInt remoteFace = 0 + nC+nV; 787 788 if (!PetscBTLookupSet(isLeaf, localFace)) { 789 localPoints[nL] = localFace; 790 remotePoints[nL].rank = neighbor; 791 remotePoints[nL].index = remoteFace; 792 } 793 } 794 #endif 795 } 796 } 797 } else if (xp > 0) { /* right */ 798 if (yp < 0) { /* bottom */ 799 if (zp < 0) { /* back */ 800 const PetscInt localVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; /* right bottom back vertex */ 801 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 802 803 if (!PetscBTLookupSet(isLeaf, localVertex)) { 804 localPoints[nL] = localVertex; 805 remotePoints[nL].rank = neighbor; 806 remotePoints[nL].index = remoteVertex; 807 ++nL; 808 } 809 } else if (zp > 0) { /* front */ 810 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; /* right bottom front vertex */ 811 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 812 813 if (!PetscBTLookupSet(isLeaf, localVertex)) { 814 localPoints[nL] = localVertex; 815 remotePoints[nL].rank = neighbor; 816 remotePoints[nL].index = remoteVertex; 817 ++nL; 818 } 819 } else { 820 nleavesCheck += nVz; 821 for (zv = 0; zv < nVz; ++zv) { 822 const PetscInt localVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; /* right bottom vertices */ 823 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 824 825 if (!PetscBTLookupSet(isLeaf, localVertex)) { 826 localPoints[nL] = localVertex; 827 remotePoints[nL].rank = neighbor; 828 remotePoints[nL].index = remoteVertex; 829 ++nL; 830 } 831 } 832 } 833 } else if (yp > 0) { /* top */ 834 if (zp < 0) { /* back */ 835 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top back vertex */ 836 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 837 838 if (!PetscBTLookupSet(isLeaf, localVertex)) { 839 localPoints[nL] = localVertex; 840 remotePoints[nL].rank = neighbor; 841 remotePoints[nL].index = remoteVertex; 842 ++nL; 843 } 844 } else if (zp > 0) { /* front */ 845 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top front vertex */ 846 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 847 848 if (!PetscBTLookupSet(isLeaf, localVertex)) { 849 localPoints[nL] = localVertex; 850 remotePoints[nL].rank = neighbor; 851 remotePoints[nL].index = remoteVertex; 852 ++nL; 853 } 854 } else { 855 for (zv = 0; zv < nVz; ++zv) { 856 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top vertices */ 857 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 858 859 if (!PetscBTLookupSet(isLeaf, localVertex)) { 860 localPoints[nL] = localVertex; 861 remotePoints[nL].rank = neighbor; 862 remotePoints[nL].index = remoteVertex; 863 ++nL; 864 } 865 } 866 } 867 } else { 868 if (zp < 0) { /* back */ 869 for (yv = 0; yv < nVy; ++yv) { 870 const PetscInt localVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; /* right back vertices */ 871 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 872 873 if (!PetscBTLookupSet(isLeaf, localVertex)) { 874 localPoints[nL] = localVertex; 875 remotePoints[nL].rank = neighbor; 876 remotePoints[nL].index = remoteVertex; 877 ++nL; 878 } 879 } 880 } else if (zp > 0) { /* front */ 881 for (yv = 0; yv < nVy; ++yv) { 882 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* right front vertices */ 883 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 884 885 if (!PetscBTLookupSet(isLeaf, localVertex)) { 886 localPoints[nL] = localVertex; 887 remotePoints[nL].rank = neighbor; 888 remotePoints[nL].index = remoteVertex; 889 ++nL; 890 } 891 } 892 } else { 893 for (zv = 0; zv < nVz; ++zv) { 894 for (yv = 0; yv < nVy; ++yv) { 895 const PetscInt localVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* right vertices */ 896 const PetscInt remoteVertex = (zv*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 897 898 if (!PetscBTLookupSet(isLeaf, localVertex)) { 899 localPoints[nL] = localVertex; 900 remotePoints[nL].rank = neighbor; 901 remotePoints[nL].index = remoteVertex; 902 ++nL; 903 } 904 } 905 } 906 #if 0 907 for (xf = 0; xf < nxF; ++xf) { 908 /* THIS IS WRONG */ 909 const PetscInt localFace = 0 + nC+nV; /* right faces */ 910 const PetscInt remoteFace = 0 + nC+nV; 911 912 if (!PetscBTLookupSet(isLeaf, localFace)) { 913 localPoints[nL] = localFace; 914 remotePoints[nL].rank = neighbor; 915 remotePoints[nL].index = remoteFace; 916 ++nL; 917 } 918 } 919 #endif 920 } 921 } 922 } else { 923 if (yp < 0) { /* bottom */ 924 if (zp < 0) { /* back */ 925 for (xv = 0; xv < nVx; ++xv) { 926 const PetscInt localVertex = ( 0*nVy + 0)*nVx + xv + nC; /* bottom back vertices */ 927 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 928 929 if (!PetscBTLookupSet(isLeaf, localVertex)) { 930 localPoints[nL] = localVertex; 931 remotePoints[nL].rank = neighbor; 932 remotePoints[nL].index = remoteVertex; 933 ++nL; 934 } 935 } 936 } else if (zp > 0) { /* front */ 937 for (xv = 0; xv < nVx; ++xv) { 938 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; /* bottom front vertices */ 939 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 940 941 if (!PetscBTLookupSet(isLeaf, localVertex)) { 942 localPoints[nL] = localVertex; 943 remotePoints[nL].rank = neighbor; 944 remotePoints[nL].index = remoteVertex; 945 ++nL; 946 } 947 } 948 } else { 949 for (zv = 0; zv < nVz; ++zv) { 950 for (xv = 0; xv < nVx; ++xv) { 951 const PetscInt localVertex = (zv*nVy + 0)*nVx + xv + nC; /* bottom vertices */ 952 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 953 954 if (!PetscBTLookupSet(isLeaf, localVertex)) { 955 localPoints[nL] = localVertex; 956 remotePoints[nL].rank = neighbor; 957 remotePoints[nL].index = remoteVertex; 958 ++nL; 959 } 960 } 961 } 962 #if 0 963 for (yf = 0; yf < nyF; ++yf) { 964 /* THIS IS WRONG */ 965 const PetscInt localFace = 0 + nC+nV; /* bottom faces */ 966 const PetscInt remoteFace = 0 + nC+nV; 967 968 if (!PetscBTLookupSet(isLeaf, localFace)) { 969 localPoints[nL] = localFace; 970 remotePoints[nL].rank = neighbor; 971 remotePoints[nL].index = remoteFace; 972 ++nL; 973 } 974 } 975 #endif 976 } 977 } else if (yp > 0) { /* top */ 978 if (zp < 0) { /* back */ 979 for (xv = 0; xv < nVx; ++xv) { 980 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; /* top back vertices */ 981 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 982 983 if (!PetscBTLookupSet(isLeaf, localVertex)) { 984 localPoints[nL] = localVertex; 985 remotePoints[nL].rank = neighbor; 986 remotePoints[nL].index = remoteVertex; 987 ++nL; 988 } 989 } 990 } else if (zp > 0) { /* front */ 991 for (xv = 0; xv < nVx; ++xv) { 992 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* top front vertices */ 993 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 994 995 if (!PetscBTLookupSet(isLeaf, localVertex)) { 996 localPoints[nL] = localVertex; 997 remotePoints[nL].rank = neighbor; 998 remotePoints[nL].index = remoteVertex; 999 ++nL; 1000 } 1001 } 1002 } else { 1003 for (zv = 0; zv < nVz; ++zv) { 1004 for (xv = 0; xv < nVx; ++xv) { 1005 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* top vertices */ 1006 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 1007 1008 if (!PetscBTLookupSet(isLeaf, localVertex)) { 1009 localPoints[nL] = localVertex; 1010 remotePoints[nL].rank = neighbor; 1011 remotePoints[nL].index = remoteVertex; 1012 ++nL; 1013 } 1014 } 1015 } 1016 #if 0 1017 for (yf = 0; yf < nyF; ++yf) { 1018 /* THIS IS WRONG */ 1019 const PetscInt localFace = 0 + nC+nV; /* top faces */ 1020 const PetscInt remoteFace = 0 + nC+nV; 1021 1022 if (!PetscBTLookupSet(isLeaf, localFace)) { 1023 localPoints[nL] = localFace; 1024 remotePoints[nL].rank = neighbor; 1025 remotePoints[nL].index = remoteFace; 1026 ++nL; 1027 } 1028 } 1029 #endif 1030 } 1031 } else { 1032 if (zp < 0) { /* back */ 1033 for (yv = 0; yv < nVy; ++yv) { 1034 for (xv = 0; xv < nVx; ++xv) { 1035 const PetscInt localVertex = ( 0*nVy + yv)*nVx + xv + nC; /* back vertices */ 1036 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 1037 1038 if (!PetscBTLookupSet(isLeaf, localVertex)) { 1039 localPoints[nL] = localVertex; 1040 remotePoints[nL].rank = neighbor; 1041 remotePoints[nL].index = remoteVertex; 1042 ++nL; 1043 } 1044 } 1045 } 1046 #if 0 1047 for (zf = 0; zf < nzF; ++zf) { 1048 /* THIS IS WRONG */ 1049 const PetscInt localFace = 0 + nC+nV; /* back faces */ 1050 const PetscInt remoteFace = 0 + nC+nV; 1051 1052 if (!PetscBTLookupSet(isLeaf, localFace)) { 1053 localPoints[nL] = localFace; 1054 remotePoints[nL].rank = neighbor; 1055 remotePoints[nL].index = remoteFace; 1056 ++nL; 1057 } 1058 } 1059 #endif 1060 } else if (zp > 0) { /* front */ 1061 for (yv = 0; yv < nVy; ++yv) { 1062 for (xv = 0; xv < nVx; ++xv) { 1063 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* front vertices */ 1064 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 1065 1066 if (!PetscBTLookupSet(isLeaf, localVertex)) { 1067 localPoints[nL] = localVertex; 1068 remotePoints[nL].rank = neighbor; 1069 remotePoints[nL].index = remoteVertex; 1070 ++nL; 1071 } 1072 } 1073 } 1074 #if 0 1075 for (zf = 0; zf < nzF; ++zf) { 1076 /* THIS IS WRONG */ 1077 const PetscInt localFace = 0 + nC+nV; /* front faces */ 1078 const PetscInt remoteFace = 0 + nC+nV; 1079 1080 if (!PetscBTLookupSet(isLeaf, localFace)) { 1081 localPoints[nL] = localFace; 1082 remotePoints[nL].rank = neighbor; 1083 remotePoints[nL].index = remoteFace; 1084 ++nL; 1085 } 1086 } 1087 #endif 1088 } else { 1089 /* Nothing is shared from the interior */ 1090 } 1091 } 1092 } 1093 } 1094 } 1095 } 1096 } 1097 ierr = PetscBTDestroy(&isLeaf);CHKERRQ(ierr); 1098 /* Remove duplication in leaf determination */ 1099 if (nleaves != nL) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "The number of leaves %d did not match the number of remote leaves %d", nleaves, nleavesCheck); 1100 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf);CHKERRQ(ierr); 1101 ierr = PetscSFSetGraph(sf, pEnd, nleaves, localPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 1102 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 1103 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 1104 *s = section; 1105 PetscFunctionReturn(0); 1106 } 1107 1108 #undef __FUNCT__ 1109 #define __FUNCT__ "DMDASetVertexCoordinates" 1110 PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu) 1111 { 1112 DM_DA *da = (DM_DA *) dm->data; 1113 Vec coordinates; 1114 PetscSection section; 1115 PetscScalar *coords; 1116 PetscReal h[3]; 1117 PetscInt dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k; 1118 PetscErrorCode ierr; 1119 1120 PetscFunctionBegin; 1121 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1122 ierr = DMDAGetInfo(dm, &dim, &M, &N, &P, 0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 1123 h[0] = (xu - xl)/M; 1124 h[1] = (yu - yl)/N; 1125 h[2] = (zu - zl)/P; 1126 ierr = DMDAGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1127 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr); 1128 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ion);CHKERRQ(ierr); 1129 ierr = PetscSectionSetNumFields(section, 1);CHKERRQ(ierr); 1130 ierr = PetscSectionSetFieldComponents(section, 0, dim);CHKERRQ(ierr); 1131 ierr = PetscSectionSetChart(section, vStart, vEnd);CHKERRQ(ierr); 1132 for (v = vStart; v < vEnd; ++v) { 1133 ierr = PetscSectionSetDof(section, v, dim);CHKERRQ(ierr); 1134 } 1135 ierr = PetscSectionSetUp(section);CHKERRQ(ierr); 1136 ierr = PetscSectionGetStorageSize(section, &size);CHKERRQ(ierr); 1137 ierr = VecCreateSeq(PETSC_COMM_SELF, size, &coordinates);CHKERRQ(ierr); 1138 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1139 for (k = 0; k < nVz; ++k) { 1140 PetscInt ind[3] = {0, 0, k + da->zs}, d, off; 1141 1142 for (j = 0; j < nVy; ++j) { 1143 ind[1] = j + da->ys; 1144 for (i = 0; i < nVx; ++i) { 1145 const PetscInt vertex = (k*nVy + j)*nVx + i + vStart; 1146 1147 ierr = PetscSectionGetOffset(section, vertex, &off);CHKERRQ(ierr); 1148 ind[0] = i + da->xs; 1149 for (d = 0; d < dim; ++d) { 1150 coords[off+d] = h[d]*ind[d]; 1151 } 1152 } 1153 } 1154 } 1155 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1156 ierr = DMSetCoordinateSection(dm, section);CHKERRQ(ierr); 1157 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1158 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr); 1159 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1160 PetscFunctionReturn(0); 1161 } 1162 PetscFunctionReturn(0); 1163 } 1164 1165 /* ------------------------------------------------------------------- */ 1166 1167 #undef __FUNCT__ 1168 #define __FUNCT__ "DMDAGetArray" 1169 /*@C 1170 DMDAGetArray - Gets a work array for a DMDA 1171 1172 Input Parameter: 1173 + da - information about my local patch 1174 - ghosted - do you want arrays for the ghosted or nonghosted patch 1175 1176 Output Parameters: 1177 . vptr - array data structured 1178 1179 Note: The vector values are NOT initialized and may have garbage in them, so you may need 1180 to zero them. 1181 1182 Level: advanced 1183 1184 .seealso: DMDARestoreArray() 1185 1186 @*/ 1187 PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr) 1188 { 1189 PetscErrorCode ierr; 1190 PetscInt j,i,xs,ys,xm,ym,zs,zm; 1191 char *iarray_start; 1192 void **iptr = (void**)vptr; 1193 DM_DA *dd = (DM_DA*)da->data; 1194 1195 PetscFunctionBegin; 1196 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1197 if (ghosted) { 1198 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1199 if (dd->arrayghostedin[i]) { 1200 *iptr = dd->arrayghostedin[i]; 1201 iarray_start = (char*)dd->startghostedin[i]; 1202 dd->arrayghostedin[i] = NULL; 1203 dd->startghostedin[i] = NULL; 1204 1205 goto done; 1206 } 1207 } 1208 xs = dd->Xs; 1209 ys = dd->Ys; 1210 zs = dd->Zs; 1211 xm = dd->Xe-dd->Xs; 1212 ym = dd->Ye-dd->Ys; 1213 zm = dd->Ze-dd->Zs; 1214 } else { 1215 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1216 if (dd->arrayin[i]) { 1217 *iptr = dd->arrayin[i]; 1218 iarray_start = (char*)dd->startin[i]; 1219 dd->arrayin[i] = NULL; 1220 dd->startin[i] = NULL; 1221 1222 goto done; 1223 } 1224 } 1225 xs = dd->xs; 1226 ys = dd->ys; 1227 zs = dd->zs; 1228 xm = dd->xe-dd->xs; 1229 ym = dd->ye-dd->ys; 1230 zm = dd->ze-dd->zs; 1231 } 1232 1233 switch (dd->dim) { 1234 case 1: { 1235 void *ptr; 1236 1237 ierr = PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1238 1239 ptr = (void*)(iarray_start - xs*sizeof(PetscScalar)); 1240 *iptr = (void*)ptr; 1241 break; 1242 } 1243 case 2: { 1244 void **ptr; 1245 1246 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1247 1248 ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*)); 1249 for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs); 1250 *iptr = (void*)ptr; 1251 break; 1252 } 1253 case 3: { 1254 void ***ptr,**bptr; 1255 1256 ierr = PetscMalloc((zm+1)*sizeof(void**)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1257 1258 ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*)); 1259 bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**)); 1260 for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys); 1261 for (i=zs; i<zs+zm; i++) { 1262 for (j=ys; j<ys+ym; j++) { 1263 ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs); 1264 } 1265 } 1266 *iptr = (void*)ptr; 1267 break; 1268 } 1269 default: 1270 SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 1271 } 1272 1273 done: 1274 /* add arrays to the checked out list */ 1275 if (ghosted) { 1276 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1277 if (!dd->arrayghostedout[i]) { 1278 dd->arrayghostedout[i] = *iptr; 1279 dd->startghostedout[i] = iarray_start; 1280 break; 1281 } 1282 } 1283 } else { 1284 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1285 if (!dd->arrayout[i]) { 1286 dd->arrayout[i] = *iptr; 1287 dd->startout[i] = iarray_start; 1288 break; 1289 } 1290 } 1291 } 1292 PetscFunctionReturn(0); 1293 } 1294 1295 #undef __FUNCT__ 1296 #define __FUNCT__ "DMDARestoreArray" 1297 /*@C 1298 DMDARestoreArray - Restores an array of derivative types for a DMDA 1299 1300 Input Parameter: 1301 + da - information about my local patch 1302 . ghosted - do you want arrays for the ghosted or nonghosted patch 1303 - vptr - array data structured to be passed to ad_FormFunctionLocal() 1304 1305 Level: advanced 1306 1307 .seealso: DMDAGetArray() 1308 1309 @*/ 1310 PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr) 1311 { 1312 PetscInt i; 1313 void **iptr = (void**)vptr,*iarray_start = 0; 1314 DM_DA *dd = (DM_DA*)da->data; 1315 1316 PetscFunctionBegin; 1317 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1318 if (ghosted) { 1319 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1320 if (dd->arrayghostedout[i] == *iptr) { 1321 iarray_start = dd->startghostedout[i]; 1322 dd->arrayghostedout[i] = NULL; 1323 dd->startghostedout[i] = NULL; 1324 break; 1325 } 1326 } 1327 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1328 if (!dd->arrayghostedin[i]) { 1329 dd->arrayghostedin[i] = *iptr; 1330 dd->startghostedin[i] = iarray_start; 1331 break; 1332 } 1333 } 1334 } else { 1335 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1336 if (dd->arrayout[i] == *iptr) { 1337 iarray_start = dd->startout[i]; 1338 dd->arrayout[i] = NULL; 1339 dd->startout[i] = NULL; 1340 break; 1341 } 1342 } 1343 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1344 if (!dd->arrayin[i]) { 1345 dd->arrayin[i] = *iptr; 1346 dd->startin[i] = iarray_start; 1347 break; 1348 } 1349 } 1350 } 1351 PetscFunctionReturn(0); 1352 } 1353 1354