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