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 h[0] = (xu - xl)/M; 1172 h[1] = (yu - yl)/N; 1173 h[2] = (zu - zl)/P; 1174 ierr = DMDAGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1175 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr); 1176 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ion);CHKERRQ(ierr); 1177 ierr = PetscSectionSetNumFields(section, 1);CHKERRQ(ierr); 1178 ierr = PetscSectionSetFieldComponents(section, 0, dim);CHKERRQ(ierr); 1179 ierr = PetscSectionSetChart(section, vStart, vEnd);CHKERRQ(ierr); 1180 for (v = vStart; v < vEnd; ++v) { 1181 ierr = PetscSectionSetDof(section, v, dim);CHKERRQ(ierr); 1182 } 1183 ierr = PetscSectionSetUp(section);CHKERRQ(ierr); 1184 ierr = PetscSectionGetStorageSize(section, &size);CHKERRQ(ierr); 1185 ierr = VecCreateSeq(PETSC_COMM_SELF, size, &coordinates);CHKERRQ(ierr); 1186 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1187 for (k = 0; k < nVz; ++k) { 1188 PetscInt ind[3], d, off; 1189 1190 ind[0] = 0; 1191 ind[1] = 0; 1192 ind[2] = k + da->zs; 1193 for (j = 0; j < nVy; ++j) { 1194 ind[1] = j + da->ys; 1195 for (i = 0; i < nVx; ++i) { 1196 const PetscInt vertex = (k*nVy + j)*nVx + i + vStart; 1197 1198 ierr = PetscSectionGetOffset(section, vertex, &off);CHKERRQ(ierr); 1199 ind[0] = i + da->xs; 1200 for (d = 0; d < dim; ++d) { 1201 coords[off+d] = h[d]*ind[d]; 1202 } 1203 } 1204 } 1205 } 1206 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1207 ierr = DMSetCoordinateSection(dm, PETSC_DETERMINE, section);CHKERRQ(ierr); 1208 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 1209 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr); 1210 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1211 PetscFunctionReturn(0); 1212 } 1213 1214 #undef __FUNCT__ 1215 #define __FUNCT__ "DMDAProjectFunctionLocal" 1216 PetscErrorCode DMDAProjectFunctionLocal(DM dm, PetscErrorCode (**funcs)(PetscInt, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 1217 { 1218 PetscDS prob; 1219 PetscFE fe; 1220 PetscDualSpace sp; 1221 PetscQuadrature q; 1222 PetscSection section; 1223 PetscScalar *values; 1224 PetscInt numFields, numComp, numPoints, dim, dimEmbed, spDim, totDim, numValues, cStart, cEnd, f, c, v, d; 1225 PetscErrorCode ierr; 1226 1227 PetscFunctionBegin; 1228 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1229 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1230 ierr = PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);CHKERRQ(ierr); 1231 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1232 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1233 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1234 ierr = DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 1235 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 1236 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1237 ierr = DMDAVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr); 1238 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim); 1239 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 1240 ierr = PetscQuadratureGetData(q, NULL, &numPoints, NULL, NULL);CHKERRQ(ierr); 1241 for (c = cStart; c < cEnd; ++c) { 1242 PetscFECellGeom geom; 1243 1244 ierr = DMDAComputeCellGeometryFEM(dm, c, q, geom.v0, geom.J, NULL, &geom.detJ);CHKERRQ(ierr); 1245 geom.dim = dim; 1246 geom.dimEmbed = dimEmbed; 1247 for (f = 0, v = 0; f < numFields; ++f) { 1248 void * const ctx = ctxs ? ctxs[f] : NULL; 1249 1250 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 1251 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 1252 ierr = PetscFEGetDualSpace(fe, &sp);CHKERRQ(ierr); 1253 ierr = PetscDualSpaceGetDimension(sp, &spDim);CHKERRQ(ierr); 1254 for (d = 0; d < spDim; ++d) { 1255 ierr = PetscDualSpaceApply(sp, d, &geom, numComp, funcs[f], ctx, &values[v]);CHKERRQ(ierr); 1256 v += numComp; 1257 } 1258 } 1259 ierr = DMDAVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr); 1260 } 1261 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 1262 PetscFunctionReturn(0); 1263 } 1264 1265 #undef __FUNCT__ 1266 #define __FUNCT__ "DMDAProjectFunction" 1267 /*@C 1268 DMDAProjectFunction - This projects the given function into the function space provided. 1269 1270 Input Parameters: 1271 + dm - The DM 1272 . funcs - The coordinate functions to evaluate 1273 . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null. 1274 - mode - The insertion mode for values 1275 1276 Output Parameter: 1277 . X - vector 1278 1279 Level: developer 1280 1281 .seealso: DMDAComputeL2Diff() 1282 @*/ 1283 PetscErrorCode DMDAProjectFunction(DM dm, PetscErrorCode (**funcs)(PetscInt, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X) 1284 { 1285 Vec localX; 1286 PetscErrorCode ierr; 1287 1288 PetscFunctionBegin; 1289 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1290 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1291 ierr = DMDAProjectFunctionLocal(dm, funcs, ctxs, mode, localX);CHKERRQ(ierr); 1292 ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr); 1293 ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr); 1294 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1295 PetscFunctionReturn(0); 1296 } 1297 1298 #undef __FUNCT__ 1299 #define __FUNCT__ "DMDAComputeL2Diff" 1300 /*@C 1301 DMDAComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h. 1302 1303 Input Parameters: 1304 + dm - The DM 1305 . funcs - The functions to evaluate for each field component 1306 . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null. 1307 - X - The coefficient vector u_h 1308 1309 Output Parameter: 1310 . diff - The diff ||u - u_h||_2 1311 1312 Level: developer 1313 1314 .seealso: DMDAProjectFunction(), DMDAComputeL2GradientDiff() 1315 @*/ 1316 PetscErrorCode DMDAComputeL2Diff(DM dm, PetscErrorCode (**funcs)(PetscInt, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 1317 { 1318 const PetscInt debug = 0; 1319 PetscDS prob; 1320 PetscFE fe; 1321 PetscQuadrature quad; 1322 PetscSection section; 1323 Vec localX; 1324 PetscScalar *funcVal; 1325 PetscReal *coords, *v0, *J, *invJ, detJ; 1326 PetscReal localDiff = 0.0; 1327 PetscInt dim, numFields, Nc, cStart, cEnd, c, field, fieldOffset, comp; 1328 PetscErrorCode ierr; 1329 1330 PetscFunctionBegin; 1331 ierr = DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 1332 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1333 ierr = PetscDSGetTotalComponents(prob, &Nc);CHKERRQ(ierr); 1334 ierr = PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);CHKERRQ(ierr); 1335 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1336 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1337 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1338 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1339 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1340 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1341 /* There are no BC values in DAs right now: ierr = DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 1342 ierr = PetscMalloc5(Nc,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1343 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1344 for (c = cStart; c < cEnd; ++c) { 1345 PetscScalar *x = NULL; 1346 PetscReal elemDiff = 0.0; 1347 1348 ierr = DMDAComputeCellGeometryFEM(dm, c, quad, v0, J, invJ, &detJ);CHKERRQ(ierr); 1349 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 1350 ierr = DMDAVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1351 1352 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 1353 void * const ctx = ctxs ? ctxs[field] : NULL; 1354 const PetscReal *quadPoints, *quadWeights; 1355 PetscReal *basis; 1356 PetscInt numQuadPoints, numBasisFuncs, numBasisComps, q, d, e, fc, f; 1357 1358 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 1359 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 1360 ierr = PetscFEGetDimension(fe, &numBasisFuncs);CHKERRQ(ierr); 1361 ierr = PetscFEGetNumComponents(fe, &numBasisComps);CHKERRQ(ierr); 1362 ierr = PetscFEGetDefaultTabulation(fe, &basis, NULL, NULL);CHKERRQ(ierr); 1363 if (debug) { 1364 char title[1024]; 1365 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 1366 ierr = DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);CHKERRQ(ierr); 1367 } 1368 for (q = 0; q < numQuadPoints; ++q) { 1369 for (d = 0; d < dim; d++) { 1370 coords[d] = v0[d]; 1371 for (e = 0; e < dim; e++) { 1372 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 1373 } 1374 } 1375 ierr = (*funcs[field])(dim, coords, numFields, funcVal, ctx);CHKERRQ(ierr); 1376 for (fc = 0; fc < numBasisComps; ++fc) { 1377 PetscScalar interpolant = 0.0; 1378 1379 for (f = 0; f < numBasisFuncs; ++f) { 1380 const PetscInt fidx = f*numBasisComps+fc; 1381 interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx]; 1382 } 1383 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);} 1384 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 1385 } 1386 } 1387 comp += numBasisComps; 1388 fieldOffset += numBasisFuncs*numBasisComps; 1389 } 1390 ierr = DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1391 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 1392 localDiff += elemDiff; 1393 } 1394 ierr = PetscFree5(funcVal,coords,v0,J,invJ);CHKERRQ(ierr); 1395 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1396 ierr = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1397 *diff = PetscSqrtReal(*diff); 1398 PetscFunctionReturn(0); 1399 } 1400 1401 #undef __FUNCT__ 1402 #define __FUNCT__ "DMDAComputeL2GradientDiff" 1403 /*@C 1404 DMDAComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h. 1405 1406 Input Parameters: 1407 + dm - The DM 1408 . funcs - The gradient functions to evaluate for each field component 1409 . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null. 1410 . X - The coefficient vector u_h 1411 - n - The vector to project along 1412 1413 Output Parameter: 1414 . diff - The diff ||(grad u - grad u_h) . n||_2 1415 1416 Level: developer 1417 1418 .seealso: DMDAProjectFunction(), DMPlexComputeL2Diff() 1419 @*/ 1420 PetscErrorCode DMDAComputeL2GradientDiff(DM dm, PetscErrorCode (**funcs)(PetscInt, const PetscReal [], const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff) 1421 { 1422 const PetscInt debug = 0; 1423 PetscDS prob; 1424 PetscFE fe; 1425 PetscQuadrature quad; 1426 PetscSection section; 1427 Vec localX; 1428 PetscScalar *funcVal, *interpolantVec; 1429 PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ; 1430 PetscReal localDiff = 0.0; 1431 PetscInt dim, numFields, Nc, cStart, cEnd, c, field, fieldOffset, comp; 1432 PetscErrorCode ierr; 1433 1434 PetscFunctionBegin; 1435 ierr = DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 1436 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1437 ierr = PetscDSGetTotalComponents(prob, &Nc);CHKERRQ(ierr); 1438 ierr = PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);CHKERRQ(ierr); 1439 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1440 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1441 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1442 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1443 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1444 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1445 /* There are no BC values in DAs right now: ierr = DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 1446 ierr = PetscMalloc7(Nc,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr); 1447 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1448 for (c = cStart; c < cEnd; ++c) { 1449 PetscScalar *x = NULL; 1450 PetscReal elemDiff = 0.0; 1451 1452 ierr = DMDAComputeCellGeometryFEM(dm, c, quad, v0, J, invJ, &detJ);CHKERRQ(ierr); 1453 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 1454 ierr = DMDAVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1455 1456 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 1457 void * const ctx = ctxs ? ctxs[field] : NULL; 1458 const PetscReal *quadPoints, *quadWeights; 1459 PetscReal *basisDer; 1460 PetscInt Nq, Nb, Nc, q, d, e, fc, f, g; 1461 1462 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 1463 ierr = PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 1464 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1465 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1466 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 1467 if (debug) { 1468 char title[1024]; 1469 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 1470 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 1471 } 1472 for (q = 0; q < Nq; ++q) { 1473 for (d = 0; d < dim; d++) { 1474 coords[d] = v0[d]; 1475 for (e = 0; e < dim; e++) { 1476 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 1477 } 1478 } 1479 ierr = (*funcs[field])(dim, coords, n, numFields, funcVal, ctx);CHKERRQ(ierr); 1480 for (fc = 0; fc < Nc; ++fc) { 1481 PetscScalar interpolant = 0.0; 1482 1483 for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0; 1484 for (f = 0; f < Nb; ++f) { 1485 const PetscInt fidx = f*Nc+fc; 1486 1487 for (d = 0; d < dim; ++d) { 1488 realSpaceDer[d] = 0.0; 1489 for (g = 0; g < dim; ++g) { 1490 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Nc+fidx)*dim+g]; 1491 } 1492 interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d]; 1493 } 1494 } 1495 for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d]; 1496 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);} 1497 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 1498 } 1499 } 1500 comp += Nc; 1501 fieldOffset += Nb*Nc; 1502 } 1503 ierr = DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1504 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 1505 localDiff += elemDiff; 1506 } 1507 ierr = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr); 1508 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1509 ierr = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1510 *diff = PetscSqrtReal(*diff); 1511 PetscFunctionReturn(0); 1512 } 1513 1514 /* ------------------------------------------------------------------- */ 1515 1516 #undef __FUNCT__ 1517 #define __FUNCT__ "DMDAGetArray" 1518 /*@C 1519 DMDAGetArray - Gets a work array for a DMDA 1520 1521 Input Parameter: 1522 + da - information about my local patch 1523 - ghosted - do you want arrays for the ghosted or nonghosted patch 1524 1525 Output Parameters: 1526 . vptr - array data structured 1527 1528 Note: The vector values are NOT initialized and may have garbage in them, so you may need 1529 to zero them. 1530 1531 Level: advanced 1532 1533 .seealso: DMDARestoreArray() 1534 1535 @*/ 1536 PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr) 1537 { 1538 PetscErrorCode ierr; 1539 PetscInt j,i,xs,ys,xm,ym,zs,zm; 1540 char *iarray_start; 1541 void **iptr = (void**)vptr; 1542 DM_DA *dd = (DM_DA*)da->data; 1543 1544 PetscFunctionBegin; 1545 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1546 if (ghosted) { 1547 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1548 if (dd->arrayghostedin[i]) { 1549 *iptr = dd->arrayghostedin[i]; 1550 iarray_start = (char*)dd->startghostedin[i]; 1551 dd->arrayghostedin[i] = NULL; 1552 dd->startghostedin[i] = NULL; 1553 1554 goto done; 1555 } 1556 } 1557 xs = dd->Xs; 1558 ys = dd->Ys; 1559 zs = dd->Zs; 1560 xm = dd->Xe-dd->Xs; 1561 ym = dd->Ye-dd->Ys; 1562 zm = dd->Ze-dd->Zs; 1563 } else { 1564 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1565 if (dd->arrayin[i]) { 1566 *iptr = dd->arrayin[i]; 1567 iarray_start = (char*)dd->startin[i]; 1568 dd->arrayin[i] = NULL; 1569 dd->startin[i] = NULL; 1570 1571 goto done; 1572 } 1573 } 1574 xs = dd->xs; 1575 ys = dd->ys; 1576 zs = dd->zs; 1577 xm = dd->xe-dd->xs; 1578 ym = dd->ye-dd->ys; 1579 zm = dd->ze-dd->zs; 1580 } 1581 1582 switch (da->dim) { 1583 case 1: { 1584 void *ptr; 1585 1586 ierr = PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1587 1588 ptr = (void*)(iarray_start - xs*sizeof(PetscScalar)); 1589 *iptr = (void*)ptr; 1590 break; 1591 } 1592 case 2: { 1593 void **ptr; 1594 1595 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1596 1597 ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*)); 1598 for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs); 1599 *iptr = (void*)ptr; 1600 break; 1601 } 1602 case 3: { 1603 void ***ptr,**bptr; 1604 1605 ierr = PetscMalloc((zm+1)*sizeof(void**)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1606 1607 ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*)); 1608 bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**)); 1609 for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys); 1610 for (i=zs; i<zs+zm; i++) { 1611 for (j=ys; j<ys+ym; j++) { 1612 ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs); 1613 } 1614 } 1615 *iptr = (void*)ptr; 1616 break; 1617 } 1618 default: 1619 SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",da->dim); 1620 } 1621 1622 done: 1623 /* add arrays to the checked out list */ 1624 if (ghosted) { 1625 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1626 if (!dd->arrayghostedout[i]) { 1627 dd->arrayghostedout[i] = *iptr; 1628 dd->startghostedout[i] = iarray_start; 1629 break; 1630 } 1631 } 1632 } else { 1633 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1634 if (!dd->arrayout[i]) { 1635 dd->arrayout[i] = *iptr; 1636 dd->startout[i] = iarray_start; 1637 break; 1638 } 1639 } 1640 } 1641 PetscFunctionReturn(0); 1642 } 1643 1644 #undef __FUNCT__ 1645 #define __FUNCT__ "DMDARestoreArray" 1646 /*@C 1647 DMDARestoreArray - Restores an array of derivative types for a DMDA 1648 1649 Input Parameter: 1650 + da - information about my local patch 1651 . ghosted - do you want arrays for the ghosted or nonghosted patch 1652 - vptr - array data structured to be passed to ad_FormFunctionLocal() 1653 1654 Level: advanced 1655 1656 .seealso: DMDAGetArray() 1657 1658 @*/ 1659 PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr) 1660 { 1661 PetscInt i; 1662 void **iptr = (void**)vptr,*iarray_start = 0; 1663 DM_DA *dd = (DM_DA*)da->data; 1664 1665 PetscFunctionBegin; 1666 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1667 if (ghosted) { 1668 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1669 if (dd->arrayghostedout[i] == *iptr) { 1670 iarray_start = dd->startghostedout[i]; 1671 dd->arrayghostedout[i] = NULL; 1672 dd->startghostedout[i] = NULL; 1673 break; 1674 } 1675 } 1676 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1677 if (!dd->arrayghostedin[i]) { 1678 dd->arrayghostedin[i] = *iptr; 1679 dd->startghostedin[i] = iarray_start; 1680 break; 1681 } 1682 } 1683 } else { 1684 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1685 if (dd->arrayout[i] == *iptr) { 1686 iarray_start = dd->startout[i]; 1687 dd->arrayout[i] = NULL; 1688 dd->startout[i] = NULL; 1689 break; 1690 } 1691 } 1692 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1693 if (!dd->arrayin[i]) { 1694 dd->arrayin[i] = *iptr; 1695 dd->startin[i] = iarray_start; 1696 break; 1697 } 1698 } 1699 } 1700 PetscFunctionReturn(0); 1701 } 1702 1703