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