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