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, const PetscInt numComp[], const PetscInt numDof[], const 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 *, void *), void **ctxs, 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, numPoints, 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 = PetscQuadratureGetData(q, NULL, &numPoints, NULL, NULL);CHKERRQ(ierr); 1247 ierr = PetscMalloc3(dim*numPoints,&v0,dim*dim*numPoints,&J,numPoints,&detJ);CHKERRQ(ierr); 1248 for (c = cStart; c < cEnd; ++c) { 1249 PetscCellGeometry geom; 1250 1251 ierr = DMDAComputeCellGeometry(dm, c, q, v0, J, NULL, detJ);CHKERRQ(ierr); 1252 geom.v0 = v0; 1253 geom.J = J; 1254 geom.detJ = detJ; 1255 for (f = 0, v = 0; f < numFields; ++f) { 1256 void * const ctx = ctxs ? ctxs[f] : NULL; 1257 ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr); 1258 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 1259 for (d = 0; d < spDim; ++d) { 1260 ierr = PetscDualSpaceApply(sp[f], d, geom, numComp, funcs[f], ctx, &values[v]);CHKERRQ(ierr); 1261 v += numComp; 1262 } 1263 } 1264 ierr = DMDAVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr); 1265 } 1266 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 1267 ierr = PetscFree3(v0,J,detJ);CHKERRQ(ierr); 1268 ierr = PetscFree(sp);CHKERRQ(ierr); 1269 PetscFunctionReturn(0); 1270 } 1271 1272 #undef __FUNCT__ 1273 #define __FUNCT__ "DMDAProjectFunction" 1274 /*@C 1275 DMDAProjectFunction - This projects the given function into the function space provided. 1276 1277 Input Parameters: 1278 + dm - The DM 1279 . fe - The PetscFE associated with the field 1280 . funcs - The coordinate functions to evaluate 1281 . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null. 1282 - mode - The insertion mode for values 1283 1284 Output Parameter: 1285 . X - vector 1286 1287 Level: developer 1288 1289 .seealso: DMDAComputeL2Diff() 1290 @*/ 1291 PetscErrorCode DMDAProjectFunction(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X) 1292 { 1293 Vec localX; 1294 PetscErrorCode ierr; 1295 1296 PetscFunctionBegin; 1297 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1298 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1299 ierr = DMDAProjectFunctionLocal(dm, fe, funcs, ctxs, mode, localX);CHKERRQ(ierr); 1300 ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr); 1301 ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr); 1302 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1303 PetscFunctionReturn(0); 1304 } 1305 1306 #undef __FUNCT__ 1307 #define __FUNCT__ "DMDAComputeL2Diff" 1308 /*@C 1309 DMDAComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h. 1310 1311 Input Parameters: 1312 + dm - The DM 1313 . fe - The PetscFE object for each field 1314 . funcs - The functions to evaluate for each field component 1315 . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null. 1316 - X - The coefficient vector u_h 1317 1318 Output Parameter: 1319 . diff - The diff ||u - u_h||_2 1320 1321 Level: developer 1322 1323 .seealso: DMDAProjectFunction() 1324 @*/ 1325 PetscErrorCode DMDAComputeL2Diff(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 1326 { 1327 const PetscInt debug = 0; 1328 PetscSection section; 1329 PetscQuadrature quad; 1330 Vec localX; 1331 PetscScalar *funcVal; 1332 PetscReal *coords, *v0, *J, *invJ, detJ; 1333 PetscReal localDiff = 0.0; 1334 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp; 1335 PetscErrorCode ierr; 1336 1337 PetscFunctionBegin; 1338 ierr = DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 1339 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1340 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1341 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1342 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1343 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1344 for (field = 0; field < numFields; ++field) { 1345 PetscInt Nc; 1346 1347 ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr); 1348 numComponents += Nc; 1349 } 1350 /* There are no BC values in DAs right now: ierr = DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 1351 ierr = PetscMalloc5(numComponents,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1352 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1353 ierr = PetscFEGetQuadrature(fe[0], &quad);CHKERRQ(ierr); 1354 for (c = cStart; c < cEnd; ++c) { 1355 PetscScalar *x = NULL; 1356 PetscReal elemDiff = 0.0; 1357 1358 ierr = DMDAComputeCellGeometry(dm, c, quad, v0, J, invJ, &detJ);CHKERRQ(ierr); 1359 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 1360 ierr = DMDAVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1361 1362 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 1363 void * const ctx = ctxs ? ctxs[field] : NULL; 1364 const PetscReal *quadPoints, *quadWeights; 1365 PetscReal *basis; 1366 PetscInt numQuadPoints, numBasisFuncs, numBasisComps, q, d, e, fc, f; 1367 1368 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 1369 ierr = PetscFEGetDimension(fe[field], &numBasisFuncs);CHKERRQ(ierr); 1370 ierr = PetscFEGetNumComponents(fe[field], &numBasisComps);CHKERRQ(ierr); 1371 ierr = PetscFEGetDefaultTabulation(fe[field], &basis, NULL, NULL);CHKERRQ(ierr); 1372 if (debug) { 1373 char title[1024]; 1374 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 1375 ierr = DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);CHKERRQ(ierr); 1376 } 1377 for (q = 0; q < numQuadPoints; ++q) { 1378 for (d = 0; d < dim; d++) { 1379 coords[d] = v0[d]; 1380 for (e = 0; e < dim; e++) { 1381 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 1382 } 1383 } 1384 (*funcs[field])(coords, funcVal, ctx); 1385 for (fc = 0; fc < numBasisComps; ++fc) { 1386 PetscScalar interpolant = 0.0; 1387 1388 for (f = 0; f < numBasisFuncs; ++f) { 1389 const PetscInt fidx = f*numBasisComps+fc; 1390 interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx]; 1391 } 1392 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);} 1393 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 1394 } 1395 } 1396 comp += numBasisComps; 1397 fieldOffset += numBasisFuncs*numBasisComps; 1398 } 1399 ierr = DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1400 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 1401 localDiff += elemDiff; 1402 } 1403 ierr = PetscFree5(funcVal,coords,v0,J,invJ);CHKERRQ(ierr); 1404 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1405 ierr = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1406 *diff = PetscSqrtReal(*diff); 1407 PetscFunctionReturn(0); 1408 } 1409 1410 /* ------------------------------------------------------------------- */ 1411 1412 #undef __FUNCT__ 1413 #define __FUNCT__ "DMDAGetArray" 1414 /*@C 1415 DMDAGetArray - Gets a work array for a DMDA 1416 1417 Input Parameter: 1418 + da - information about my local patch 1419 - ghosted - do you want arrays for the ghosted or nonghosted patch 1420 1421 Output Parameters: 1422 . vptr - array data structured 1423 1424 Note: The vector values are NOT initialized and may have garbage in them, so you may need 1425 to zero them. 1426 1427 Level: advanced 1428 1429 .seealso: DMDARestoreArray() 1430 1431 @*/ 1432 PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr) 1433 { 1434 PetscErrorCode ierr; 1435 PetscInt j,i,xs,ys,xm,ym,zs,zm; 1436 char *iarray_start; 1437 void **iptr = (void**)vptr; 1438 DM_DA *dd = (DM_DA*)da->data; 1439 1440 PetscFunctionBegin; 1441 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1442 if (ghosted) { 1443 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1444 if (dd->arrayghostedin[i]) { 1445 *iptr = dd->arrayghostedin[i]; 1446 iarray_start = (char*)dd->startghostedin[i]; 1447 dd->arrayghostedin[i] = NULL; 1448 dd->startghostedin[i] = NULL; 1449 1450 goto done; 1451 } 1452 } 1453 xs = dd->Xs; 1454 ys = dd->Ys; 1455 zs = dd->Zs; 1456 xm = dd->Xe-dd->Xs; 1457 ym = dd->Ye-dd->Ys; 1458 zm = dd->Ze-dd->Zs; 1459 } else { 1460 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1461 if (dd->arrayin[i]) { 1462 *iptr = dd->arrayin[i]; 1463 iarray_start = (char*)dd->startin[i]; 1464 dd->arrayin[i] = NULL; 1465 dd->startin[i] = NULL; 1466 1467 goto done; 1468 } 1469 } 1470 xs = dd->xs; 1471 ys = dd->ys; 1472 zs = dd->zs; 1473 xm = dd->xe-dd->xs; 1474 ym = dd->ye-dd->ys; 1475 zm = dd->ze-dd->zs; 1476 } 1477 1478 switch (dd->dim) { 1479 case 1: { 1480 void *ptr; 1481 1482 ierr = PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1483 1484 ptr = (void*)(iarray_start - xs*sizeof(PetscScalar)); 1485 *iptr = (void*)ptr; 1486 break; 1487 } 1488 case 2: { 1489 void **ptr; 1490 1491 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1492 1493 ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*)); 1494 for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs); 1495 *iptr = (void*)ptr; 1496 break; 1497 } 1498 case 3: { 1499 void ***ptr,**bptr; 1500 1501 ierr = PetscMalloc((zm+1)*sizeof(void**)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1502 1503 ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*)); 1504 bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**)); 1505 for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys); 1506 for (i=zs; i<zs+zm; i++) { 1507 for (j=ys; j<ys+ym; j++) { 1508 ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs); 1509 } 1510 } 1511 *iptr = (void*)ptr; 1512 break; 1513 } 1514 default: 1515 SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 1516 } 1517 1518 done: 1519 /* add arrays to the checked out list */ 1520 if (ghosted) { 1521 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1522 if (!dd->arrayghostedout[i]) { 1523 dd->arrayghostedout[i] = *iptr; 1524 dd->startghostedout[i] = iarray_start; 1525 break; 1526 } 1527 } 1528 } else { 1529 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1530 if (!dd->arrayout[i]) { 1531 dd->arrayout[i] = *iptr; 1532 dd->startout[i] = iarray_start; 1533 break; 1534 } 1535 } 1536 } 1537 PetscFunctionReturn(0); 1538 } 1539 1540 #undef __FUNCT__ 1541 #define __FUNCT__ "DMDARestoreArray" 1542 /*@C 1543 DMDARestoreArray - Restores an array of derivative types for a DMDA 1544 1545 Input Parameter: 1546 + da - information about my local patch 1547 . ghosted - do you want arrays for the ghosted or nonghosted patch 1548 - vptr - array data structured to be passed to ad_FormFunctionLocal() 1549 1550 Level: advanced 1551 1552 .seealso: DMDAGetArray() 1553 1554 @*/ 1555 PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr) 1556 { 1557 PetscInt i; 1558 void **iptr = (void**)vptr,*iarray_start = 0; 1559 DM_DA *dd = (DM_DA*)da->data; 1560 1561 PetscFunctionBegin; 1562 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1563 if (ghosted) { 1564 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1565 if (dd->arrayghostedout[i] == *iptr) { 1566 iarray_start = dd->startghostedout[i]; 1567 dd->arrayghostedout[i] = NULL; 1568 dd->startghostedout[i] = NULL; 1569 break; 1570 } 1571 } 1572 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1573 if (!dd->arrayghostedin[i]) { 1574 dd->arrayghostedin[i] = *iptr; 1575 dd->startghostedin[i] = iarray_start; 1576 break; 1577 } 1578 } 1579 } else { 1580 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1581 if (dd->arrayout[i] == *iptr) { 1582 iarray_start = dd->startout[i]; 1583 dd->arrayout[i] = NULL; 1584 dd->startout[i] = NULL; 1585 break; 1586 } 1587 } 1588 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1589 if (!dd->arrayin[i]) { 1590 dd->arrayin[i] = *iptr; 1591 dd->startin[i] = iarray_start; 1592 break; 1593 } 1594 } 1595 } 1596 PetscFunctionReturn(0); 1597 } 1598 1599