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