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