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 = PetscMalloc2(nleaves,&localPoints,nleaves,&remotePoints);CHKERRQ(ierr); 713 for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) { 714 for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) { 715 for (xn = 0; xn < 3; ++xn) { 716 const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0; 717 const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn]; 718 PetscInt xv, yv, zv; 719 720 if (neighbor >= 0 && neighbor < rank) { 721 if (xp < 0) { /* left */ 722 if (yp < 0) { /* bottom */ 723 if (zp < 0) { /* back */ 724 const PetscInt localVertex = ( 0*nVy + 0)*nVx + 0 + nC; /* left bottom back vertex */ 725 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 726 727 if (!PetscBTLookupSet(isLeaf, localVertex)) { 728 localPoints[nL] = localVertex; 729 remotePoints[nL].rank = neighbor; 730 remotePoints[nL].index = remoteVertex; 731 ++nL; 732 } 733 } else if (zp > 0) { /* front */ 734 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; /* left bottom front vertex */ 735 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 736 737 if (!PetscBTLookupSet(isLeaf, localVertex)) { 738 localPoints[nL] = localVertex; 739 remotePoints[nL].rank = neighbor; 740 remotePoints[nL].index = remoteVertex; 741 ++nL; 742 } 743 } else { 744 for (zv = 0; zv < nVz; ++zv) { 745 const PetscInt localVertex = (zv*nVy + 0)*nVx + 0 + nC; /* left bottom vertices */ 746 const PetscInt remoteVertex = (zv*nVy + nVy-1)*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 } 755 } 756 } else if (yp > 0) { /* top */ 757 if (zp < 0) { /* back */ 758 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; /* left top back vertex */ 759 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 760 761 if (!PetscBTLookupSet(isLeaf, localVertex)) { 762 localPoints[nL] = localVertex; 763 remotePoints[nL].rank = neighbor; 764 remotePoints[nL].index = remoteVertex; 765 ++nL; 766 } 767 } else if (zp > 0) { /* front */ 768 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; /* left top front vertex */ 769 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 770 771 if (!PetscBTLookupSet(isLeaf, localVertex)) { 772 localPoints[nL] = localVertex; 773 remotePoints[nL].rank = neighbor; 774 remotePoints[nL].index = remoteVertex; 775 ++nL; 776 } 777 } else { 778 for (zv = 0; zv < nVz; ++zv) { 779 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; /* left top vertices */ 780 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 781 782 if (!PetscBTLookupSet(isLeaf, localVertex)) { 783 localPoints[nL] = localVertex; 784 remotePoints[nL].rank = neighbor; 785 remotePoints[nL].index = remoteVertex; 786 ++nL; 787 } 788 } 789 } 790 } else { 791 if (zp < 0) { /* back */ 792 for (yv = 0; yv < nVy; ++yv) { 793 const PetscInt localVertex = ( 0*nVy + yv)*nVx + 0 + nC; /* left back vertices */ 794 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 795 796 if (!PetscBTLookupSet(isLeaf, localVertex)) { 797 localPoints[nL] = localVertex; 798 remotePoints[nL].rank = neighbor; 799 remotePoints[nL].index = remoteVertex; 800 ++nL; 801 } 802 } 803 } else if (zp > 0) { /* front */ 804 for (yv = 0; yv < nVy; ++yv) { 805 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; /* left front vertices */ 806 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 807 808 if (!PetscBTLookupSet(isLeaf, localVertex)) { 809 localPoints[nL] = localVertex; 810 remotePoints[nL].rank = neighbor; 811 remotePoints[nL].index = remoteVertex; 812 ++nL; 813 } 814 } 815 } else { 816 for (zv = 0; zv < nVz; ++zv) { 817 for (yv = 0; yv < nVy; ++yv) { 818 const PetscInt localVertex = (zv*nVy + yv)*nVx + 0 + nC; /* left vertices */ 819 const PetscInt remoteVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 820 821 if (!PetscBTLookupSet(isLeaf, localVertex)) { 822 localPoints[nL] = localVertex; 823 remotePoints[nL].rank = neighbor; 824 remotePoints[nL].index = remoteVertex; 825 ++nL; 826 } 827 } 828 } 829 #if 0 830 for (xf = 0; xf < nxF; ++xf) { 831 /* THIS IS WRONG */ 832 const PetscInt localFace = 0 + nC+nV; /* left faces */ 833 const PetscInt remoteFace = 0 + nC+nV; 834 835 if (!PetscBTLookupSet(isLeaf, localFace)) { 836 localPoints[nL] = localFace; 837 remotePoints[nL].rank = neighbor; 838 remotePoints[nL].index = remoteFace; 839 } 840 } 841 #endif 842 } 843 } 844 } else if (xp > 0) { /* right */ 845 if (yp < 0) { /* bottom */ 846 if (zp < 0) { /* back */ 847 const PetscInt localVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; /* right bottom back vertex */ 848 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 849 850 if (!PetscBTLookupSet(isLeaf, localVertex)) { 851 localPoints[nL] = localVertex; 852 remotePoints[nL].rank = neighbor; 853 remotePoints[nL].index = remoteVertex; 854 ++nL; 855 } 856 } else if (zp > 0) { /* front */ 857 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; /* right bottom front vertex */ 858 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 859 860 if (!PetscBTLookupSet(isLeaf, localVertex)) { 861 localPoints[nL] = localVertex; 862 remotePoints[nL].rank = neighbor; 863 remotePoints[nL].index = remoteVertex; 864 ++nL; 865 } 866 } else { 867 nleavesCheck += nVz; 868 for (zv = 0; zv < nVz; ++zv) { 869 const PetscInt localVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; /* right bottom vertices */ 870 const PetscInt remoteVertex = (zv*nVy + nVy-1)*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 } 879 } 880 } else if (yp > 0) { /* top */ 881 if (zp < 0) { /* back */ 882 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top back vertex */ 883 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 884 885 if (!PetscBTLookupSet(isLeaf, localVertex)) { 886 localPoints[nL] = localVertex; 887 remotePoints[nL].rank = neighbor; 888 remotePoints[nL].index = remoteVertex; 889 ++nL; 890 } 891 } else if (zp > 0) { /* front */ 892 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top front vertex */ 893 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 894 895 if (!PetscBTLookupSet(isLeaf, localVertex)) { 896 localPoints[nL] = localVertex; 897 remotePoints[nL].rank = neighbor; 898 remotePoints[nL].index = remoteVertex; 899 ++nL; 900 } 901 } else { 902 for (zv = 0; zv < nVz; ++zv) { 903 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top vertices */ 904 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 905 906 if (!PetscBTLookupSet(isLeaf, localVertex)) { 907 localPoints[nL] = localVertex; 908 remotePoints[nL].rank = neighbor; 909 remotePoints[nL].index = remoteVertex; 910 ++nL; 911 } 912 } 913 } 914 } else { 915 if (zp < 0) { /* back */ 916 for (yv = 0; yv < nVy; ++yv) { 917 const PetscInt localVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; /* right back vertices */ 918 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 919 920 if (!PetscBTLookupSet(isLeaf, localVertex)) { 921 localPoints[nL] = localVertex; 922 remotePoints[nL].rank = neighbor; 923 remotePoints[nL].index = remoteVertex; 924 ++nL; 925 } 926 } 927 } else if (zp > 0) { /* front */ 928 for (yv = 0; yv < nVy; ++yv) { 929 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* right front vertices */ 930 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 931 932 if (!PetscBTLookupSet(isLeaf, localVertex)) { 933 localPoints[nL] = localVertex; 934 remotePoints[nL].rank = neighbor; 935 remotePoints[nL].index = remoteVertex; 936 ++nL; 937 } 938 } 939 } else { 940 for (zv = 0; zv < nVz; ++zv) { 941 for (yv = 0; yv < nVy; ++yv) { 942 const PetscInt localVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* right vertices */ 943 const PetscInt remoteVertex = (zv*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 944 945 if (!PetscBTLookupSet(isLeaf, localVertex)) { 946 localPoints[nL] = localVertex; 947 remotePoints[nL].rank = neighbor; 948 remotePoints[nL].index = remoteVertex; 949 ++nL; 950 } 951 } 952 } 953 #if 0 954 for (xf = 0; xf < nxF; ++xf) { 955 /* THIS IS WRONG */ 956 const PetscInt localFace = 0 + nC+nV; /* right faces */ 957 const PetscInt remoteFace = 0 + nC+nV; 958 959 if (!PetscBTLookupSet(isLeaf, localFace)) { 960 localPoints[nL] = localFace; 961 remotePoints[nL].rank = neighbor; 962 remotePoints[nL].index = remoteFace; 963 ++nL; 964 } 965 } 966 #endif 967 } 968 } 969 } else { 970 if (yp < 0) { /* bottom */ 971 if (zp < 0) { /* back */ 972 for (xv = 0; xv < nVx; ++xv) { 973 const PetscInt localVertex = ( 0*nVy + 0)*nVx + xv + nC; /* bottom back vertices */ 974 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 975 976 if (!PetscBTLookupSet(isLeaf, localVertex)) { 977 localPoints[nL] = localVertex; 978 remotePoints[nL].rank = neighbor; 979 remotePoints[nL].index = remoteVertex; 980 ++nL; 981 } 982 } 983 } else if (zp > 0) { /* front */ 984 for (xv = 0; xv < nVx; ++xv) { 985 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; /* bottom front vertices */ 986 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 987 988 if (!PetscBTLookupSet(isLeaf, localVertex)) { 989 localPoints[nL] = localVertex; 990 remotePoints[nL].rank = neighbor; 991 remotePoints[nL].index = remoteVertex; 992 ++nL; 993 } 994 } 995 } else { 996 for (zv = 0; zv < nVz; ++zv) { 997 for (xv = 0; xv < nVx; ++xv) { 998 const PetscInt localVertex = (zv*nVy + 0)*nVx + xv + nC; /* bottom vertices */ 999 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 1000 1001 if (!PetscBTLookupSet(isLeaf, localVertex)) { 1002 localPoints[nL] = localVertex; 1003 remotePoints[nL].rank = neighbor; 1004 remotePoints[nL].index = remoteVertex; 1005 ++nL; 1006 } 1007 } 1008 } 1009 #if 0 1010 for (yf = 0; yf < nyF; ++yf) { 1011 /* THIS IS WRONG */ 1012 const PetscInt localFace = 0 + nC+nV; /* bottom faces */ 1013 const PetscInt remoteFace = 0 + nC+nV; 1014 1015 if (!PetscBTLookupSet(isLeaf, localFace)) { 1016 localPoints[nL] = localFace; 1017 remotePoints[nL].rank = neighbor; 1018 remotePoints[nL].index = remoteFace; 1019 ++nL; 1020 } 1021 } 1022 #endif 1023 } 1024 } else if (yp > 0) { /* top */ 1025 if (zp < 0) { /* back */ 1026 for (xv = 0; xv < nVx; ++xv) { 1027 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; /* top back vertices */ 1028 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 1029 1030 if (!PetscBTLookupSet(isLeaf, localVertex)) { 1031 localPoints[nL] = localVertex; 1032 remotePoints[nL].rank = neighbor; 1033 remotePoints[nL].index = remoteVertex; 1034 ++nL; 1035 } 1036 } 1037 } else if (zp > 0) { /* front */ 1038 for (xv = 0; xv < nVx; ++xv) { 1039 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* top front vertices */ 1040 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 1041 1042 if (!PetscBTLookupSet(isLeaf, localVertex)) { 1043 localPoints[nL] = localVertex; 1044 remotePoints[nL].rank = neighbor; 1045 remotePoints[nL].index = remoteVertex; 1046 ++nL; 1047 } 1048 } 1049 } else { 1050 for (zv = 0; zv < nVz; ++zv) { 1051 for (xv = 0; xv < nVx; ++xv) { 1052 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* top vertices */ 1053 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 1054 1055 if (!PetscBTLookupSet(isLeaf, localVertex)) { 1056 localPoints[nL] = localVertex; 1057 remotePoints[nL].rank = neighbor; 1058 remotePoints[nL].index = remoteVertex; 1059 ++nL; 1060 } 1061 } 1062 } 1063 #if 0 1064 for (yf = 0; yf < nyF; ++yf) { 1065 /* THIS IS WRONG */ 1066 const PetscInt localFace = 0 + nC+nV; /* top faces */ 1067 const PetscInt remoteFace = 0 + nC+nV; 1068 1069 if (!PetscBTLookupSet(isLeaf, localFace)) { 1070 localPoints[nL] = localFace; 1071 remotePoints[nL].rank = neighbor; 1072 remotePoints[nL].index = remoteFace; 1073 ++nL; 1074 } 1075 } 1076 #endif 1077 } 1078 } else { 1079 if (zp < 0) { /* back */ 1080 for (yv = 0; yv < nVy; ++yv) { 1081 for (xv = 0; xv < nVx; ++xv) { 1082 const PetscInt localVertex = ( 0*nVy + yv)*nVx + xv + nC; /* back vertices */ 1083 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 1084 1085 if (!PetscBTLookupSet(isLeaf, localVertex)) { 1086 localPoints[nL] = localVertex; 1087 remotePoints[nL].rank = neighbor; 1088 remotePoints[nL].index = remoteVertex; 1089 ++nL; 1090 } 1091 } 1092 } 1093 #if 0 1094 for (zf = 0; zf < nzF; ++zf) { 1095 /* THIS IS WRONG */ 1096 const PetscInt localFace = 0 + nC+nV; /* back faces */ 1097 const PetscInt remoteFace = 0 + nC+nV; 1098 1099 if (!PetscBTLookupSet(isLeaf, localFace)) { 1100 localPoints[nL] = localFace; 1101 remotePoints[nL].rank = neighbor; 1102 remotePoints[nL].index = remoteFace; 1103 ++nL; 1104 } 1105 } 1106 #endif 1107 } else if (zp > 0) { /* front */ 1108 for (yv = 0; yv < nVy; ++yv) { 1109 for (xv = 0; xv < nVx; ++xv) { 1110 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* front vertices */ 1111 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 1112 1113 if (!PetscBTLookupSet(isLeaf, localVertex)) { 1114 localPoints[nL] = localVertex; 1115 remotePoints[nL].rank = neighbor; 1116 remotePoints[nL].index = remoteVertex; 1117 ++nL; 1118 } 1119 } 1120 } 1121 #if 0 1122 for (zf = 0; zf < nzF; ++zf) { 1123 /* THIS IS WRONG */ 1124 const PetscInt localFace = 0 + nC+nV; /* front faces */ 1125 const PetscInt remoteFace = 0 + nC+nV; 1126 1127 if (!PetscBTLookupSet(isLeaf, localFace)) { 1128 localPoints[nL] = localFace; 1129 remotePoints[nL].rank = neighbor; 1130 remotePoints[nL].index = remoteFace; 1131 ++nL; 1132 } 1133 } 1134 #endif 1135 } else { 1136 /* Nothing is shared from the interior */ 1137 } 1138 } 1139 } 1140 } 1141 } 1142 } 1143 } 1144 ierr = PetscBTDestroy(&isLeaf);CHKERRQ(ierr); 1145 /* Remove duplication in leaf determination */ 1146 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); 1147 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf);CHKERRQ(ierr); 1148 ierr = PetscSFSetGraph(sf, pEnd, nleaves, localPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 1149 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 1150 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 1151 *s = section; 1152 PetscFunctionReturn(0); 1153 } 1154 1155 #undef __FUNCT__ 1156 #define __FUNCT__ "DMDASetVertexCoordinates" 1157 PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu) 1158 { 1159 DM_DA *da = (DM_DA *) dm->data; 1160 Vec coordinates; 1161 PetscSection section; 1162 PetscScalar *coords; 1163 PetscReal h[3]; 1164 PetscInt dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k; 1165 PetscErrorCode ierr; 1166 1167 PetscFunctionBegin; 1168 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1169 ierr = DMDAGetInfo(dm, &dim, &M, &N, &P, 0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 1170 h[0] = (xu - xl)/M; 1171 h[1] = (yu - yl)/N; 1172 h[2] = (zu - zl)/P; 1173 ierr = DMDAGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1174 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr); 1175 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ion);CHKERRQ(ierr); 1176 ierr = PetscSectionSetNumFields(section, 1);CHKERRQ(ierr); 1177 ierr = PetscSectionSetFieldComponents(section, 0, dim);CHKERRQ(ierr); 1178 ierr = PetscSectionSetChart(section, vStart, vEnd);CHKERRQ(ierr); 1179 for (v = vStart; v < vEnd; ++v) { 1180 ierr = PetscSectionSetDof(section, v, dim);CHKERRQ(ierr); 1181 } 1182 ierr = PetscSectionSetUp(section);CHKERRQ(ierr); 1183 ierr = PetscSectionGetStorageSize(section, &size);CHKERRQ(ierr); 1184 ierr = VecCreateSeq(PETSC_COMM_SELF, size, &coordinates);CHKERRQ(ierr); 1185 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1186 for (k = 0; k < nVz; ++k) { 1187 PetscInt ind[3], d, off; 1188 1189 ind[0] = 0; 1190 ind[1] = 0; 1191 ind[2] = k + da->zs; 1192 for (j = 0; j < nVy; ++j) { 1193 ind[1] = j + da->ys; 1194 for (i = 0; i < nVx; ++i) { 1195 const PetscInt vertex = (k*nVy + j)*nVx + i + vStart; 1196 1197 ierr = PetscSectionGetOffset(section, vertex, &off);CHKERRQ(ierr); 1198 ind[0] = i + da->xs; 1199 for (d = 0; d < dim; ++d) { 1200 coords[off+d] = h[d]*ind[d]; 1201 } 1202 } 1203 } 1204 } 1205 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1206 ierr = DMSetCoordinateSection(dm, PETSC_DETERMINE, section);CHKERRQ(ierr); 1207 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1208 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr); 1209 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1210 PetscFunctionReturn(0); 1211 } 1212 1213 #undef __FUNCT__ 1214 #define __FUNCT__ "DMDAProjectFunctionLocal" 1215 PetscErrorCode DMDAProjectFunctionLocal(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 1216 { 1217 PetscDS prob; 1218 PetscFE fe; 1219 PetscDualSpace sp; 1220 PetscQuadrature q; 1221 PetscSection section; 1222 PetscScalar *values; 1223 PetscInt numFields, numComp, numPoints, dim, spDim, totDim, numValues, cStart, cEnd, f, c, v, d; 1224 PetscErrorCode ierr; 1225 1226 PetscFunctionBegin; 1227 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1228 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1229 ierr = PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);CHKERRQ(ierr); 1230 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1231 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1232 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1233 ierr = DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 1234 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1235 ierr = DMDAVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr); 1236 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim); 1237 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 1238 ierr = PetscQuadratureGetData(q, NULL, &numPoints, NULL, NULL);CHKERRQ(ierr); 1239 for (c = cStart; c < cEnd; ++c) { 1240 PetscFECellGeom geom; 1241 1242 ierr = DMDAComputeCellGeometryFEM(dm, c, q, geom.v0, geom.J, NULL, &geom.detJ);CHKERRQ(ierr); 1243 for (f = 0, v = 0; f < numFields; ++f) { 1244 void * const ctx = ctxs ? ctxs[f] : NULL; 1245 1246 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 1247 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 1248 ierr = PetscFEGetDualSpace(fe, &sp);CHKERRQ(ierr); 1249 ierr = PetscDualSpaceGetDimension(sp, &spDim);CHKERRQ(ierr); 1250 for (d = 0; d < spDim; ++d) { 1251 ierr = PetscDualSpaceApply(sp, d, &geom, numComp, funcs[f], ctx, &values[v]);CHKERRQ(ierr); 1252 v += numComp; 1253 } 1254 } 1255 ierr = DMDAVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr); 1256 } 1257 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 1258 PetscFunctionReturn(0); 1259 } 1260 1261 #undef __FUNCT__ 1262 #define __FUNCT__ "DMDAProjectFunction" 1263 /*@C 1264 DMDAProjectFunction - This projects the given function into the function space provided. 1265 1266 Input Parameters: 1267 + dm - The DM 1268 . funcs - The coordinate functions to evaluate 1269 . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null. 1270 - mode - The insertion mode for values 1271 1272 Output Parameter: 1273 . X - vector 1274 1275 Level: developer 1276 1277 .seealso: DMDAComputeL2Diff() 1278 @*/ 1279 PetscErrorCode DMDAProjectFunction(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X) 1280 { 1281 Vec localX; 1282 PetscErrorCode ierr; 1283 1284 PetscFunctionBegin; 1285 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1286 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1287 ierr = DMDAProjectFunctionLocal(dm, funcs, ctxs, mode, localX);CHKERRQ(ierr); 1288 ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr); 1289 ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr); 1290 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1291 PetscFunctionReturn(0); 1292 } 1293 1294 #undef __FUNCT__ 1295 #define __FUNCT__ "DMDAComputeL2Diff" 1296 /*@C 1297 DMDAComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h. 1298 1299 Input Parameters: 1300 + dm - The DM 1301 . funcs - The functions to evaluate for each field component 1302 . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null. 1303 - X - The coefficient vector u_h 1304 1305 Output Parameter: 1306 . diff - The diff ||u - u_h||_2 1307 1308 Level: developer 1309 1310 .seealso: DMDAProjectFunction(), DMDAComputeL2GradientDiff() 1311 @*/ 1312 PetscErrorCode DMDAComputeL2Diff(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 1313 { 1314 const PetscInt debug = 0; 1315 PetscDS prob; 1316 PetscFE fe; 1317 PetscQuadrature quad; 1318 PetscSection section; 1319 Vec localX; 1320 PetscScalar *funcVal; 1321 PetscReal *coords, *v0, *J, *invJ, detJ; 1322 PetscReal localDiff = 0.0; 1323 PetscInt dim, numFields, Nc, cStart, cEnd, c, field, fieldOffset, comp; 1324 PetscErrorCode ierr; 1325 1326 PetscFunctionBegin; 1327 ierr = DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 1328 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1329 ierr = PetscDSGetTotalComponents(prob, &Nc);CHKERRQ(ierr); 1330 ierr = PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);CHKERRQ(ierr); 1331 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1332 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1333 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1334 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1335 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1336 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1337 /* There are no BC values in DAs right now: ierr = DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 1338 ierr = PetscMalloc5(Nc,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1339 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1340 for (c = cStart; c < cEnd; ++c) { 1341 PetscScalar *x = NULL; 1342 PetscReal elemDiff = 0.0; 1343 1344 ierr = DMDAComputeCellGeometryFEM(dm, c, quad, v0, J, invJ, &detJ);CHKERRQ(ierr); 1345 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 1346 ierr = DMDAVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1347 1348 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 1349 void * const ctx = ctxs ? ctxs[field] : NULL; 1350 const PetscReal *quadPoints, *quadWeights; 1351 PetscReal *basis; 1352 PetscInt numQuadPoints, numBasisFuncs, numBasisComps, q, d, e, fc, f; 1353 1354 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 1355 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 1356 ierr = PetscFEGetDimension(fe, &numBasisFuncs);CHKERRQ(ierr); 1357 ierr = PetscFEGetNumComponents(fe, &numBasisComps);CHKERRQ(ierr); 1358 ierr = PetscFEGetDefaultTabulation(fe, &basis, NULL, NULL);CHKERRQ(ierr); 1359 if (debug) { 1360 char title[1024]; 1361 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 1362 ierr = DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);CHKERRQ(ierr); 1363 } 1364 for (q = 0; q < numQuadPoints; ++q) { 1365 for (d = 0; d < dim; d++) { 1366 coords[d] = v0[d]; 1367 for (e = 0; e < dim; e++) { 1368 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 1369 } 1370 } 1371 (*funcs[field])(coords, funcVal, ctx); 1372 for (fc = 0; fc < numBasisComps; ++fc) { 1373 PetscScalar interpolant = 0.0; 1374 1375 for (f = 0; f < numBasisFuncs; ++f) { 1376 const PetscInt fidx = f*numBasisComps+fc; 1377 interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx]; 1378 } 1379 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);} 1380 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 1381 } 1382 } 1383 comp += numBasisComps; 1384 fieldOffset += numBasisFuncs*numBasisComps; 1385 } 1386 ierr = DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1387 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 1388 localDiff += elemDiff; 1389 } 1390 ierr = PetscFree5(funcVal,coords,v0,J,invJ);CHKERRQ(ierr); 1391 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1392 ierr = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1393 *diff = PetscSqrtReal(*diff); 1394 PetscFunctionReturn(0); 1395 } 1396 1397 #undef __FUNCT__ 1398 #define __FUNCT__ "DMDAComputeL2GradientDiff" 1399 /*@C 1400 DMDAComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h. 1401 1402 Input Parameters: 1403 + dm - The DM 1404 . funcs - The gradient functions to evaluate for each field component 1405 . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null. 1406 . X - The coefficient vector u_h 1407 - n - The vector to project along 1408 1409 Output Parameter: 1410 . diff - The diff ||(grad u - grad u_h) . n||_2 1411 1412 Level: developer 1413 1414 .seealso: DMDAProjectFunction(), DMPlexComputeL2Diff() 1415 @*/ 1416 PetscErrorCode DMDAComputeL2GradientDiff(DM dm, void (**funcs)(const PetscReal [], const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff) 1417 { 1418 const PetscInt debug = 0; 1419 PetscDS prob; 1420 PetscFE fe; 1421 PetscQuadrature quad; 1422 PetscSection section; 1423 Vec localX; 1424 PetscScalar *funcVal, *interpolantVec; 1425 PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ; 1426 PetscReal localDiff = 0.0; 1427 PetscInt dim, numFields, Nc, cStart, cEnd, c, field, fieldOffset, comp; 1428 PetscErrorCode ierr; 1429 1430 PetscFunctionBegin; 1431 ierr = DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 1432 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1433 ierr = PetscDSGetTotalComponents(prob, &Nc);CHKERRQ(ierr); 1434 ierr = PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);CHKERRQ(ierr); 1435 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1436 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1437 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1438 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1439 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1440 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1441 /* There are no BC values in DAs right now: ierr = DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 1442 ierr = PetscMalloc7(Nc,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr); 1443 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1444 for (c = cStart; c < cEnd; ++c) { 1445 PetscScalar *x = NULL; 1446 PetscReal elemDiff = 0.0; 1447 1448 ierr = DMDAComputeCellGeometryFEM(dm, c, quad, v0, J, invJ, &detJ);CHKERRQ(ierr); 1449 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 1450 ierr = DMDAVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1451 1452 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 1453 void * const ctx = ctxs ? ctxs[field] : NULL; 1454 const PetscReal *quadPoints, *quadWeights; 1455 PetscReal *basisDer; 1456 PetscInt Nq, Nb, Nc, q, d, e, fc, f, g; 1457 1458 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 1459 ierr = PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 1460 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1461 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1462 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 1463 if (debug) { 1464 char title[1024]; 1465 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 1466 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 1467 } 1468 for (q = 0; q < Nq; ++q) { 1469 for (d = 0; d < dim; d++) { 1470 coords[d] = v0[d]; 1471 for (e = 0; e < dim; e++) { 1472 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 1473 } 1474 } 1475 (*funcs[field])(coords, n, funcVal, ctx); 1476 for (fc = 0; fc < Nc; ++fc) { 1477 PetscScalar interpolant = 0.0; 1478 1479 for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0; 1480 for (f = 0; f < Nb; ++f) { 1481 const PetscInt fidx = f*Nc+fc; 1482 1483 for (d = 0; d < dim; ++d) { 1484 realSpaceDer[d] = 0.0; 1485 for (g = 0; g < dim; ++g) { 1486 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Nc+fidx)*dim+g]; 1487 } 1488 interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d]; 1489 } 1490 } 1491 for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d]; 1492 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);} 1493 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 1494 } 1495 } 1496 comp += Nc; 1497 fieldOffset += Nb*Nc; 1498 } 1499 ierr = DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1500 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 1501 localDiff += elemDiff; 1502 } 1503 ierr = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr); 1504 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1505 ierr = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1506 *diff = PetscSqrtReal(*diff); 1507 PetscFunctionReturn(0); 1508 } 1509 1510 /* ------------------------------------------------------------------- */ 1511 1512 #undef __FUNCT__ 1513 #define __FUNCT__ "DMDAGetArray" 1514 /*@C 1515 DMDAGetArray - Gets a work array for a DMDA 1516 1517 Input Parameter: 1518 + da - information about my local patch 1519 - ghosted - do you want arrays for the ghosted or nonghosted patch 1520 1521 Output Parameters: 1522 . vptr - array data structured 1523 1524 Note: The vector values are NOT initialized and may have garbage in them, so you may need 1525 to zero them. 1526 1527 Level: advanced 1528 1529 .seealso: DMDARestoreArray() 1530 1531 @*/ 1532 PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr) 1533 { 1534 PetscErrorCode ierr; 1535 PetscInt j,i,xs,ys,xm,ym,zs,zm; 1536 char *iarray_start; 1537 void **iptr = (void**)vptr; 1538 DM_DA *dd = (DM_DA*)da->data; 1539 1540 PetscFunctionBegin; 1541 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1542 if (ghosted) { 1543 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1544 if (dd->arrayghostedin[i]) { 1545 *iptr = dd->arrayghostedin[i]; 1546 iarray_start = (char*)dd->startghostedin[i]; 1547 dd->arrayghostedin[i] = NULL; 1548 dd->startghostedin[i] = NULL; 1549 1550 goto done; 1551 } 1552 } 1553 xs = dd->Xs; 1554 ys = dd->Ys; 1555 zs = dd->Zs; 1556 xm = dd->Xe-dd->Xs; 1557 ym = dd->Ye-dd->Ys; 1558 zm = dd->Ze-dd->Zs; 1559 } else { 1560 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1561 if (dd->arrayin[i]) { 1562 *iptr = dd->arrayin[i]; 1563 iarray_start = (char*)dd->startin[i]; 1564 dd->arrayin[i] = NULL; 1565 dd->startin[i] = NULL; 1566 1567 goto done; 1568 } 1569 } 1570 xs = dd->xs; 1571 ys = dd->ys; 1572 zs = dd->zs; 1573 xm = dd->xe-dd->xs; 1574 ym = dd->ye-dd->ys; 1575 zm = dd->ze-dd->zs; 1576 } 1577 1578 switch (da->dim) { 1579 case 1: { 1580 void *ptr; 1581 1582 ierr = PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1583 1584 ptr = (void*)(iarray_start - xs*sizeof(PetscScalar)); 1585 *iptr = (void*)ptr; 1586 break; 1587 } 1588 case 2: { 1589 void **ptr; 1590 1591 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1592 1593 ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*)); 1594 for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs); 1595 *iptr = (void*)ptr; 1596 break; 1597 } 1598 case 3: { 1599 void ***ptr,**bptr; 1600 1601 ierr = PetscMalloc((zm+1)*sizeof(void**)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1602 1603 ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*)); 1604 bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**)); 1605 for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys); 1606 for (i=zs; i<zs+zm; i++) { 1607 for (j=ys; j<ys+ym; j++) { 1608 ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs); 1609 } 1610 } 1611 *iptr = (void*)ptr; 1612 break; 1613 } 1614 default: 1615 SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",da->dim); 1616 } 1617 1618 done: 1619 /* add arrays to the checked out list */ 1620 if (ghosted) { 1621 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1622 if (!dd->arrayghostedout[i]) { 1623 dd->arrayghostedout[i] = *iptr; 1624 dd->startghostedout[i] = iarray_start; 1625 break; 1626 } 1627 } 1628 } else { 1629 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1630 if (!dd->arrayout[i]) { 1631 dd->arrayout[i] = *iptr; 1632 dd->startout[i] = iarray_start; 1633 break; 1634 } 1635 } 1636 } 1637 PetscFunctionReturn(0); 1638 } 1639 1640 #undef __FUNCT__ 1641 #define __FUNCT__ "DMDARestoreArray" 1642 /*@C 1643 DMDARestoreArray - Restores an array of derivative types for a DMDA 1644 1645 Input Parameter: 1646 + da - information about my local patch 1647 . ghosted - do you want arrays for the ghosted or nonghosted patch 1648 - vptr - array data structured to be passed to ad_FormFunctionLocal() 1649 1650 Level: advanced 1651 1652 .seealso: DMDAGetArray() 1653 1654 @*/ 1655 PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr) 1656 { 1657 PetscInt i; 1658 void **iptr = (void**)vptr,*iarray_start = 0; 1659 DM_DA *dd = (DM_DA*)da->data; 1660 1661 PetscFunctionBegin; 1662 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1663 if (ghosted) { 1664 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1665 if (dd->arrayghostedout[i] == *iptr) { 1666 iarray_start = dd->startghostedout[i]; 1667 dd->arrayghostedout[i] = NULL; 1668 dd->startghostedout[i] = NULL; 1669 break; 1670 } 1671 } 1672 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1673 if (!dd->arrayghostedin[i]) { 1674 dd->arrayghostedin[i] = *iptr; 1675 dd->startghostedin[i] = iarray_start; 1676 break; 1677 } 1678 } 1679 } else { 1680 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1681 if (dd->arrayout[i] == *iptr) { 1682 iarray_start = dd->startout[i]; 1683 dd->arrayout[i] = NULL; 1684 dd->startout[i] = NULL; 1685 break; 1686 } 1687 } 1688 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1689 if (!dd->arrayin[i]) { 1690 dd->arrayin[i] = *iptr; 1691 dd->startin[i] = iarray_start; 1692 break; 1693 } 1694 } 1695 } 1696 PetscFunctionReturn(0); 1697 } 1698 1699