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