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