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 <petscsf.h> 8 #include <petscfe.h> 9 10 /* 11 This allows the DMDA vectors to properly tell MATLAB their dimensions 12 */ 13 #if defined(PETSC_HAVE_MATLAB_ENGINE) 14 #include <engine.h> /* MATLAB include file */ 15 #include <mex.h> /* MATLAB include file */ 16 #undef __FUNCT__ 17 #define __FUNCT__ "VecMatlabEnginePut_DA2d" 18 static PetscErrorCode VecMatlabEnginePut_DA2d(PetscObject obj,void *mengine) 19 { 20 PetscErrorCode ierr; 21 PetscInt n,m; 22 Vec vec = (Vec)obj; 23 PetscScalar *array; 24 mxArray *mat; 25 DM da; 26 27 PetscFunctionBegin; 28 ierr = VecGetDM(vec, &da);CHKERRQ(ierr); 29 if (!da) SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_ARG_WRONGSTATE,"Vector not associated with a DMDA"); 30 ierr = DMDAGetGhostCorners(da,0,0,0,&m,&n,0);CHKERRQ(ierr); 31 32 ierr = VecGetArray(vec,&array);CHKERRQ(ierr); 33 #if !defined(PETSC_USE_COMPLEX) 34 mat = mxCreateDoubleMatrix(m,n,mxREAL); 35 #else 36 mat = mxCreateDoubleMatrix(m,n,mxCOMPLEX); 37 #endif 38 ierr = PetscMemcpy(mxGetPr(mat),array,n*m*sizeof(PetscScalar));CHKERRQ(ierr); 39 ierr = PetscObjectName(obj);CHKERRQ(ierr); 40 engPutVariable((Engine*)mengine,obj->name,mat); 41 42 ierr = VecRestoreArray(vec,&array);CHKERRQ(ierr); 43 PetscFunctionReturn(0); 44 } 45 #endif 46 47 48 #undef __FUNCT__ 49 #define __FUNCT__ "DMCreateLocalVector_DA" 50 PetscErrorCode DMCreateLocalVector_DA(DM da,Vec *g) 51 { 52 PetscErrorCode ierr; 53 DM_DA *dd = (DM_DA*)da->data; 54 55 PetscFunctionBegin; 56 PetscValidHeaderSpecific(da,DM_CLASSID,1); 57 PetscValidPointer(g,2); 58 if (da->defaultSection) { 59 ierr = DMCreateLocalVector_Section_Private(da,g);CHKERRQ(ierr); 60 } else { 61 ierr = VecCreate(PETSC_COMM_SELF,g);CHKERRQ(ierr); 62 ierr = VecSetSizes(*g,dd->nlocal,PETSC_DETERMINE);CHKERRQ(ierr); 63 ierr = VecSetBlockSize(*g,dd->w);CHKERRQ(ierr); 64 ierr = VecSetType(*g,da->vectype);CHKERRQ(ierr); 65 ierr = VecSetDM(*g, da);CHKERRQ(ierr); 66 #if defined(PETSC_HAVE_MATLAB_ENGINE) 67 if (dd->w == 1 && dd->dim == 2) { 68 ierr = PetscObjectComposeFunction((PetscObject)*g,"PetscMatlabEnginePut_C",VecMatlabEnginePut_DA2d);CHKERRQ(ierr); 69 } 70 #endif 71 } 72 PetscFunctionReturn(0); 73 } 74 75 #undef __FUNCT__ 76 #define __FUNCT__ "DMDAGetNumCells" 77 PetscErrorCode DMDAGetNumCells(DM dm, PetscInt *numCellsX, PetscInt *numCellsY, PetscInt *numCellsZ, PetscInt *numCells) 78 { 79 DM_DA *da = (DM_DA*) dm->data; 80 const PetscInt dim = da->dim; 81 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 82 const PetscInt nC = (mx)*(dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1); 83 84 PetscFunctionBegin; 85 if (numCellsX) { 86 PetscValidIntPointer(numCellsX,2); 87 *numCellsX = mx; 88 } 89 if (numCellsY) { 90 PetscValidIntPointer(numCellsX,3); 91 *numCellsY = my; 92 } 93 if (numCellsZ) { 94 PetscValidIntPointer(numCellsX,4); 95 *numCellsZ = mz; 96 } 97 if (numCells) { 98 PetscValidIntPointer(numCells,5); 99 *numCells = nC; 100 } 101 PetscFunctionReturn(0); 102 } 103 104 #undef __FUNCT__ 105 #define __FUNCT__ "DMDAGetNumVertices" 106 PetscErrorCode DMDAGetNumVertices(DM dm, PetscInt *numVerticesX, PetscInt *numVerticesY, PetscInt *numVerticesZ, PetscInt *numVertices) 107 { 108 DM_DA *da = (DM_DA*) dm->data; 109 const PetscInt dim = da->dim; 110 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 111 const PetscInt nVx = mx+1; 112 const PetscInt nVy = dim > 1 ? (my+1) : 1; 113 const PetscInt nVz = dim > 2 ? (mz+1) : 1; 114 const PetscInt nV = nVx*nVy*nVz; 115 116 PetscFunctionBegin; 117 if (numVerticesX) { 118 PetscValidIntPointer(numVerticesX,2); 119 *numVerticesX = nVx; 120 } 121 if (numVerticesY) { 122 PetscValidIntPointer(numVerticesY,3); 123 *numVerticesY = nVy; 124 } 125 if (numVerticesZ) { 126 PetscValidIntPointer(numVerticesZ,4); 127 *numVerticesZ = nVz; 128 } 129 if (numVertices) { 130 PetscValidIntPointer(numVertices,5); 131 *numVertices = nV; 132 } 133 PetscFunctionReturn(0); 134 } 135 136 #undef __FUNCT__ 137 #define __FUNCT__ "DMDAGetNumFaces" 138 PetscErrorCode DMDAGetNumFaces(DM dm, PetscInt *numXFacesX, PetscInt *numXFaces, PetscInt *numYFacesY, PetscInt *numYFaces, PetscInt *numZFacesZ, PetscInt *numZFaces) 139 { 140 DM_DA *da = (DM_DA*) dm->data; 141 const PetscInt dim = da->dim; 142 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 143 const PetscInt nxF = (dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1); 144 const PetscInt nXF = (mx+1)*nxF; 145 const PetscInt nyF = mx*(dim > 2 ? mz : 1); 146 const PetscInt nYF = dim > 1 ? (my+1)*nyF : 0; 147 const PetscInt nzF = mx*(dim > 1 ? my : 0); 148 const PetscInt nZF = dim > 2 ? (mz+1)*nzF : 0; 149 150 PetscFunctionBegin; 151 if (numXFacesX) { 152 PetscValidIntPointer(numXFacesX,2); 153 *numXFacesX = nxF; 154 } 155 if (numXFaces) { 156 PetscValidIntPointer(numXFaces,3); 157 *numXFaces = nXF; 158 } 159 if (numYFacesY) { 160 PetscValidIntPointer(numYFacesY,4); 161 *numYFacesY = nyF; 162 } 163 if (numYFaces) { 164 PetscValidIntPointer(numYFaces,5); 165 *numYFaces = nYF; 166 } 167 if (numZFacesZ) { 168 PetscValidIntPointer(numZFacesZ,6); 169 *numZFacesZ = nzF; 170 } 171 if (numZFaces) { 172 PetscValidIntPointer(numZFaces,7); 173 *numZFaces = nZF; 174 } 175 PetscFunctionReturn(0); 176 } 177 178 #undef __FUNCT__ 179 #define __FUNCT__ "DMDAGetHeightStratum" 180 PetscErrorCode DMDAGetHeightStratum(DM dm, PetscInt height, PetscInt *pStart, PetscInt *pEnd) 181 { 182 DM_DA *da = (DM_DA*) dm->data; 183 const PetscInt dim = da->dim; 184 PetscInt nC, nV, nXF, nYF, nZF; 185 PetscErrorCode ierr; 186 187 PetscFunctionBegin; 188 if (pStart) PetscValidIntPointer(pStart,3); 189 if (pEnd) PetscValidIntPointer(pEnd,4); 190 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); 191 ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr); 192 ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr); 193 if (height == 0) { 194 /* Cells */ 195 if (pStart) *pStart = 0; 196 if (pEnd) *pEnd = nC; 197 } else if (height == 1) { 198 /* Faces */ 199 if (pStart) *pStart = nC+nV; 200 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 201 } else if (height == dim) { 202 /* Vertices */ 203 if (pStart) *pStart = nC; 204 if (pEnd) *pEnd = nC+nV; 205 } else if (height < 0) { 206 /* All points */ 207 if (pStart) *pStart = 0; 208 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 209 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of height %d in the DA", height); 210 PetscFunctionReturn(0); 211 } 212 213 #undef __FUNCT__ 214 #define __FUNCT__ "DMDAGetDepthStratum" 215 PetscErrorCode DMDAGetDepthStratum(DM dm, PetscInt depth, PetscInt *pStart, PetscInt *pEnd) 216 { 217 DM_DA *da = (DM_DA*) dm->data; 218 const PetscInt dim = da->dim; 219 PetscInt nC, nV, nXF, nYF, nZF; 220 PetscErrorCode ierr; 221 222 PetscFunctionBegin; 223 if (pStart) PetscValidIntPointer(pStart,3); 224 if (pEnd) PetscValidIntPointer(pEnd,4); 225 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); 226 ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr); 227 ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr); 228 if (depth == dim) { 229 /* Cells */ 230 if (pStart) *pStart = 0; 231 if (pEnd) *pEnd = nC; 232 } else if (depth == dim-1) { 233 /* Faces */ 234 if (pStart) *pStart = nC+nV; 235 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 236 } else if (depth == 0) { 237 /* Vertices */ 238 if (pStart) *pStart = nC; 239 if (pEnd) *pEnd = nC+nV; 240 } else if (depth < 0) { 241 /* All points */ 242 if (pStart) *pStart = 0; 243 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 244 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of depth %d in the DA", depth); 245 PetscFunctionReturn(0); 246 } 247 248 #undef __FUNCT__ 249 #define __FUNCT__ "DMDAGetConeSize" 250 PetscErrorCode DMDAGetConeSize(DM dm, PetscInt p, PetscInt *coneSize) 251 { 252 DM_DA *da = (DM_DA*) dm->data; 253 const PetscInt dim = da->dim; 254 PetscInt nC, nV, nXF, nYF, nZF; 255 PetscErrorCode ierr; 256 257 PetscFunctionBegin; 258 *coneSize = 0; 259 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); 260 ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr); 261 ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr); 262 switch (dim) { 263 case 2: 264 if (p >= 0) { 265 if (p < nC) { 266 *coneSize = 4; 267 } else if (p < nC+nV) { 268 *coneSize = 0; 269 } else if (p < nC+nV+nXF+nYF+nZF) { 270 *coneSize = 2; 271 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF); 272 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p); 273 break; 274 case 3: 275 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D"); 276 break; 277 } 278 PetscFunctionReturn(0); 279 } 280 281 #undef __FUNCT__ 282 #define __FUNCT__ "DMDAGetCone" 283 PetscErrorCode DMDAGetCone(DM dm, PetscInt p, PetscInt *cone[]) 284 { 285 DM_DA *da = (DM_DA*) dm->data; 286 const PetscInt dim = da->dim; 287 PetscInt nCx, nCy, nCz, nC, nVx, nVy, nVz, nV, nxF, nyF, nzF, nXF, nYF, nZF; 288 PetscErrorCode ierr; 289 290 PetscFunctionBegin; 291 if (!cone) {ierr = DMGetWorkArray(dm, 6, PETSC_INT, cone);CHKERRQ(ierr);} 292 ierr = DMDAGetNumCells(dm, &nCx, &nCy, &nCz, &nC);CHKERRQ(ierr); 293 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr); 294 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 295 switch (dim) { 296 case 2: 297 if (p >= 0) { 298 if (p < nC) { 299 const PetscInt cy = p / nCx; 300 const PetscInt cx = p % nCx; 301 302 (*cone)[0] = cy*nxF + cx + nC+nV; 303 (*cone)[1] = cx*nyF + cy + nyF + nC+nV+nXF; 304 (*cone)[2] = cy*nxF + cx + nxF + nC+nV; 305 (*cone)[3] = cx*nyF + cy + nC+nV+nXF; 306 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do cell cones"); 307 } else if (p < nC+nV) { 308 } else if (p < nC+nV+nXF) { 309 const PetscInt fy = (p - nC+nV) / nxF; 310 const PetscInt fx = (p - nC+nV) % nxF; 311 312 (*cone)[0] = fy*nVx + fx + nC; 313 (*cone)[1] = fy*nVx + fx + 1 + nC; 314 } else if (p < nC+nV+nXF+nYF) { 315 const PetscInt fx = (p - nC+nV+nXF) / nyF; 316 const PetscInt fy = (p - nC+nV+nXF) % nyF; 317 318 (*cone)[0] = fy*nVx + fx + nC; 319 (*cone)[1] = fy*nVx + fx + nVx + nC; 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__ "DMDARestoreCone" 332 PetscErrorCode DMDARestoreCone(DM dm, PetscInt p, PetscInt *cone[]) 333 { 334 PetscErrorCode ierr; 335 336 PetscFunctionBegin; 337 ierr = DMGetWorkArray(dm, 6, PETSC_INT, cone);CHKERRQ(ierr); 338 PetscFunctionReturn(0); 339 } 340 341 #undef __FUNCT__ 342 #define __FUNCT__ "DMDACreateSection" 343 /*@C 344 DMDACreateSection - Create a PetscSection inside the DMDA that describes data layout. This allows multiple fields with 345 different numbers of dofs on vertices, cells, and faces in each direction. 346 347 Input Parameters: 348 + dm- The DMDA 349 . numFields - The number of fields 350 . numComp - The number of components in each field 351 . numDof - The number of dofs per dimension for each field 352 . numFaceDof - The number of dofs per face for each field and direction, or NULL 353 354 Level: developer 355 356 Note: 357 The default DMDA numbering is as follows: 358 359 - Cells: [0, nC) 360 - Vertices: [nC, nC+nV) 361 - X-Faces: [nC+nV, nC+nV+nXF) normal is +- x-dir 362 - Y-Faces: [nC+nV+nXF, nC+nV+nXF+nYF) normal is +- y-dir 363 - Z-Faces: [nC+nV+nXF+nYF, nC+nV+nXF+nYF+nZF) normal is +- z-dir 364 365 We interpret the default DMDA partition as a cell partition, and the data assignment as a cell assignment. 366 @*/ 367 PetscErrorCode DMDACreateSection(DM dm, PetscInt numComp[], PetscInt numDof[], PetscInt numFaceDof[], PetscSection *s) 368 { 369 DM_DA *da = (DM_DA*) dm->data; 370 PetscSection section; 371 const PetscInt dim = da->dim; 372 PetscInt numFields, numVertexTotDof = 0, numCellTotDof = 0, numFaceTotDof[3] = {0, 0, 0}; 373 PetscSF sf; 374 PetscMPIInt rank; 375 const PetscMPIInt *neighbors; 376 PetscInt *localPoints; 377 PetscSFNode *remotePoints; 378 PetscInt nleaves = 0, nleavesCheck = 0, nL = 0; 379 PetscInt nC, nVx, nVy, nVz, nV, nxF, nXF, nyF, nYF, nzF, nZF; 380 PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart, zfEnd; 381 PetscInt f, v, c, xf, yf, zf, xn, yn, zn; 382 PetscErrorCode ierr; 383 384 PetscFunctionBegin; 385 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 386 PetscValidPointer(s, 4); 387 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 388 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); 389 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr); 390 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 391 ierr = DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);CHKERRQ(ierr); 392 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 393 ierr = DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 394 ierr = DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);CHKERRQ(ierr); 395 xfStart = vEnd; xfEnd = xfStart+nXF; 396 yfStart = xfEnd; yfEnd = yfStart+nYF; 397 zfStart = yfEnd; zfEnd = zfStart+nZF; 398 if (zfEnd != fEnd) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid face end %d, should be %d", zfEnd, fEnd); 399 /* Create local section */ 400 ierr = DMDAGetInfo(dm, 0,0,0,0,0,0,0, &numFields, 0,0,0,0,0);CHKERRQ(ierr); 401 for (f = 0; f < numFields; ++f) { 402 numVertexTotDof += numDof[f*(dim+1)+0]; 403 numCellTotDof += numDof[f*(dim+1)+dim]; 404 numFaceTotDof[0] += dim > 0 ? (numFaceDof ? numFaceDof[f*dim+0] : numDof[f*(dim+1)+dim-1]) : 0; 405 numFaceTotDof[1] += dim > 1 ? (numFaceDof ? numFaceDof[f*dim+1] : numDof[f*(dim+1)+dim-1]) : 0; 406 numFaceTotDof[2] += dim > 2 ? (numFaceDof ? numFaceDof[f*dim+2] : numDof[f*(dim+1)+dim-1]) : 0; 407 } 408 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion);CHKERRQ(ierr); 409 if (numFields > 0) { 410 ierr = PetscSectionSetNumFields(section, numFields);CHKERRQ(ierr); 411 for (f = 0; f < numFields; ++f) { 412 const char *name; 413 414 ierr = DMDAGetFieldName(dm, f, &name);CHKERRQ(ierr); 415 ierr = PetscSectionSetFieldName(section, f, name ? name : "Field");CHKERRQ(ierr); 416 if (numComp) { 417 ierr = PetscSectionSetFieldComponents(section, f, numComp[f]);CHKERRQ(ierr); 418 } 419 } 420 } 421 ierr = PetscSectionSetChart(section, pStart, pEnd);CHKERRQ(ierr); 422 for (v = vStart; v < vEnd; ++v) { 423 for (f = 0; f < numFields; ++f) { 424 ierr = PetscSectionSetFieldDof(section, v, f, numDof[f*(dim+1)+0]);CHKERRQ(ierr); 425 } 426 ierr = PetscSectionSetDof(section, v, numVertexTotDof);CHKERRQ(ierr); 427 } 428 for (xf = xfStart; xf < xfEnd; ++xf) { 429 for (f = 0; f < numFields; ++f) { 430 ierr = PetscSectionSetFieldDof(section, xf, f, numFaceDof ? numFaceDof[f*dim+0] : numDof[f*(dim+1)+dim-1]);CHKERRQ(ierr); 431 } 432 ierr = PetscSectionSetDof(section, xf, numFaceTotDof[0]);CHKERRQ(ierr); 433 } 434 for (yf = yfStart; yf < yfEnd; ++yf) { 435 for (f = 0; f < numFields; ++f) { 436 ierr = PetscSectionSetFieldDof(section, yf, f, numFaceDof ? numFaceDof[f*dim+1] : numDof[f*(dim+1)+dim-1]);CHKERRQ(ierr); 437 } 438 ierr = PetscSectionSetDof(section, yf, numFaceTotDof[1]);CHKERRQ(ierr); 439 } 440 for (zf = zfStart; zf < zfEnd; ++zf) { 441 for (f = 0; f < numFields; ++f) { 442 ierr = PetscSectionSetFieldDof(section, zf, f, numFaceDof ? numFaceDof[f*dim+2] : numDof[f*(dim+1)+dim-1]);CHKERRQ(ierr); 443 } 444 ierr = PetscSectionSetDof(section, zf, numFaceTotDof[2]);CHKERRQ(ierr); 445 } 446 for (c = cStart; c < cEnd; ++c) { 447 for (f = 0; f < numFields; ++f) { 448 ierr = PetscSectionSetFieldDof(section, c, f, numDof[f*(dim+1)+dim]);CHKERRQ(ierr); 449 } 450 ierr = PetscSectionSetDof(section, c, numCellTotDof);CHKERRQ(ierr); 451 } 452 ierr = PetscSectionSetUp(section);CHKERRQ(ierr); 453 /* Create mesh point SF */ 454 ierr = DMDAGetNeighbors(dm, &neighbors);CHKERRQ(ierr); 455 for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) { 456 for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) { 457 for (xn = 0; xn < 3; ++xn) { 458 const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0; 459 const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn]; 460 461 if (neighbor >= 0 && neighbor < rank) { 462 nleaves += (!xp ? nVx : 1) * (!yp ? nVy : 1) * (!zp ? nVz : 1); /* vertices */ 463 if (xp && !yp && !zp) { 464 nleaves += nxF; /* x faces */ 465 } else if (yp && !zp && !xp) { 466 nleaves += nyF; /* y faces */ 467 } else if (zp && !xp && !yp) { 468 nleaves += nzF; /* z faces */ 469 } 470 } 471 } 472 } 473 } 474 ierr = PetscMalloc2(nleaves,&localPoints,nleaves,&remotePoints);CHKERRQ(ierr); 475 for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) { 476 for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) { 477 for (xn = 0; xn < 3; ++xn) { 478 const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0; 479 const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn]; 480 PetscInt xv, yv, zv; 481 482 if (neighbor >= 0 && neighbor < rank) { 483 if (xp < 0) { /* left */ 484 if (yp < 0) { /* bottom */ 485 if (zp < 0) { /* back */ 486 const PetscInt localVertex = ( 0*nVy + 0)*nVx + 0 + nC; 487 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 488 nleavesCheck += 1; /* left bottom back vertex */ 489 490 localPoints[nL] = localVertex; 491 remotePoints[nL].rank = neighbor; 492 remotePoints[nL].index = remoteVertex; 493 ++nL; 494 } else if (zp > 0) { /* front */ 495 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; 496 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 497 nleavesCheck += 1; /* left bottom front vertex */ 498 499 localPoints[nL] = localVertex; 500 remotePoints[nL].rank = neighbor; 501 remotePoints[nL].index = remoteVertex; 502 ++nL; 503 } else { 504 nleavesCheck += nVz; /* left bottom vertices */ 505 for (zv = 0; zv < nVz; ++zv, ++nL) { 506 const PetscInt localVertex = (zv*nVy + 0)*nVx + 0 + nC; 507 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 508 509 localPoints[nL] = localVertex; 510 remotePoints[nL].rank = neighbor; 511 remotePoints[nL].index = remoteVertex; 512 } 513 } 514 } else if (yp > 0) { /* top */ 515 if (zp < 0) { /* back */ 516 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; 517 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 518 nleavesCheck += 1; /* left top back vertex */ 519 520 localPoints[nL] = localVertex; 521 remotePoints[nL].rank = neighbor; 522 remotePoints[nL].index = remoteVertex; 523 ++nL; 524 } else if (zp > 0) { /* front */ 525 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; 526 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 527 nleavesCheck += 1; /* left top front vertex */ 528 529 localPoints[nL] = localVertex; 530 remotePoints[nL].rank = neighbor; 531 remotePoints[nL].index = remoteVertex; 532 ++nL; 533 } else { 534 nleavesCheck += nVz; /* left top vertices */ 535 for (zv = 0; zv < nVz; ++zv, ++nL) { 536 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; 537 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 538 539 localPoints[nL] = localVertex; 540 remotePoints[nL].rank = neighbor; 541 remotePoints[nL].index = remoteVertex; 542 } 543 } 544 } else { 545 if (zp < 0) { /* back */ 546 nleavesCheck += nVy; /* left back vertices */ 547 for (yv = 0; yv < nVy; ++yv, ++nL) { 548 const PetscInt localVertex = ( 0*nVy + yv)*nVx + 0 + nC; 549 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 550 551 localPoints[nL] = localVertex; 552 remotePoints[nL].rank = neighbor; 553 remotePoints[nL].index = remoteVertex; 554 } 555 } else if (zp > 0) { /* front */ 556 nleavesCheck += nVy; /* left front vertices */ 557 for (yv = 0; yv < nVy; ++yv, ++nL) { 558 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; 559 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 560 561 localPoints[nL] = localVertex; 562 remotePoints[nL].rank = neighbor; 563 remotePoints[nL].index = remoteVertex; 564 } 565 } else { 566 nleavesCheck += nVy*nVz; /* left vertices */ 567 for (zv = 0; zv < nVz; ++zv) { 568 for (yv = 0; yv < nVy; ++yv, ++nL) { 569 const PetscInt localVertex = (zv*nVy + yv)*nVx + 0 + nC; 570 const PetscInt remoteVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 571 572 localPoints[nL] = localVertex; 573 remotePoints[nL].rank = neighbor; 574 remotePoints[nL].index = remoteVertex; 575 } 576 } 577 nleavesCheck += nxF; /* left faces */ 578 for (xf = 0; xf < nxF; ++xf, ++nL) { 579 /* THIS IS WRONG */ 580 const PetscInt localFace = 0 + nC+nV; 581 const PetscInt remoteFace = 0 + nC+nV; 582 583 localPoints[nL] = localFace; 584 remotePoints[nL].rank = neighbor; 585 remotePoints[nL].index = remoteFace; 586 } 587 } 588 } 589 } else if (xp > 0) { /* right */ 590 if (yp < 0) { /* bottom */ 591 if (zp < 0) { /* back */ 592 const PetscInt localVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; 593 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 594 nleavesCheck += 1; /* right bottom back vertex */ 595 596 localPoints[nL] = localVertex; 597 remotePoints[nL].rank = neighbor; 598 remotePoints[nL].index = remoteVertex; 599 ++nL; 600 } else if (zp > 0) { /* front */ 601 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; 602 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 603 nleavesCheck += 1; /* right bottom front vertex */ 604 605 localPoints[nL] = localVertex; 606 remotePoints[nL].rank = neighbor; 607 remotePoints[nL].index = remoteVertex; 608 ++nL; 609 } else { 610 nleavesCheck += nVz; /* right bottom vertices */ 611 for (zv = 0; zv < nVz; ++zv, ++nL) { 612 const PetscInt localVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; 613 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 614 615 localPoints[nL] = localVertex; 616 remotePoints[nL].rank = neighbor; 617 remotePoints[nL].index = remoteVertex; 618 } 619 } 620 } else if (yp > 0) { /* top */ 621 if (zp < 0) { /* back */ 622 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; 623 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 624 nleavesCheck += 1; /* right top back vertex */ 625 626 localPoints[nL] = localVertex; 627 remotePoints[nL].rank = neighbor; 628 remotePoints[nL].index = remoteVertex; 629 ++nL; 630 } else if (zp > 0) { /* front */ 631 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; 632 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 633 nleavesCheck += 1; /* right top front vertex */ 634 635 localPoints[nL] = localVertex; 636 remotePoints[nL].rank = neighbor; 637 remotePoints[nL].index = remoteVertex; 638 ++nL; 639 } else { 640 nleavesCheck += nVz; /* right top vertices */ 641 for (zv = 0; zv < nVz; ++zv, ++nL) { 642 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; 643 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 644 645 localPoints[nL] = localVertex; 646 remotePoints[nL].rank = neighbor; 647 remotePoints[nL].index = remoteVertex; 648 } 649 } 650 } else { 651 if (zp < 0) { /* back */ 652 nleavesCheck += nVy; /* right back vertices */ 653 for (yv = 0; yv < nVy; ++yv, ++nL) { 654 const PetscInt localVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; 655 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 656 657 localPoints[nL] = localVertex; 658 remotePoints[nL].rank = neighbor; 659 remotePoints[nL].index = remoteVertex; 660 } 661 } else if (zp > 0) { /* front */ 662 nleavesCheck += nVy; /* right front vertices */ 663 for (yv = 0; yv < nVy; ++yv, ++nL) { 664 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; 665 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 666 667 localPoints[nL] = localVertex; 668 remotePoints[nL].rank = neighbor; 669 remotePoints[nL].index = remoteVertex; 670 } 671 } else { 672 nleavesCheck += nVy*nVz; /* right vertices */ 673 for (zv = 0; zv < nVz; ++zv) { 674 for (yv = 0; yv < nVy; ++yv, ++nL) { 675 const PetscInt localVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; 676 const PetscInt remoteVertex = (zv*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 677 678 localPoints[nL] = localVertex; 679 remotePoints[nL].rank = neighbor; 680 remotePoints[nL].index = remoteVertex; 681 } 682 } 683 nleavesCheck += nxF; /* right faces */ 684 for (xf = 0; xf < nxF; ++xf, ++nL) { 685 /* THIS IS WRONG */ 686 const PetscInt localFace = 0 + nC+nV; 687 const PetscInt remoteFace = 0 + nC+nV; 688 689 localPoints[nL] = localFace; 690 remotePoints[nL].rank = neighbor; 691 remotePoints[nL].index = remoteFace; 692 } 693 } 694 } 695 } else { 696 if (yp < 0) { /* bottom */ 697 if (zp < 0) { /* back */ 698 nleavesCheck += nVx; /* bottom back vertices */ 699 for (xv = 0; xv < nVx; ++xv, ++nL) { 700 const PetscInt localVertex = ( 0*nVy + 0)*nVx + xv + nC; 701 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 702 703 localPoints[nL] = localVertex; 704 remotePoints[nL].rank = neighbor; 705 remotePoints[nL].index = remoteVertex; 706 } 707 } else if (zp > 0) { /* front */ 708 nleavesCheck += nVx; /* bottom front vertices */ 709 for (xv = 0; xv < nVx; ++xv, ++nL) { 710 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; 711 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 712 713 localPoints[nL] = localVertex; 714 remotePoints[nL].rank = neighbor; 715 remotePoints[nL].index = remoteVertex; 716 } 717 } else { 718 nleavesCheck += nVx*nVz; /* bottom vertices */ 719 for (zv = 0; zv < nVz; ++zv) { 720 for (xv = 0; xv < nVx; ++xv, ++nL) { 721 const PetscInt localVertex = (zv*nVy + 0)*nVx + xv + nC; 722 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 723 724 localPoints[nL] = localVertex; 725 remotePoints[nL].rank = neighbor; 726 remotePoints[nL].index = remoteVertex; 727 } 728 } 729 nleavesCheck += nyF; /* bottom faces */ 730 for (yf = 0; yf < nyF; ++yf, ++nL) { 731 /* THIS IS WRONG */ 732 const PetscInt localFace = 0 + nC+nV; 733 const PetscInt remoteFace = 0 + nC+nV; 734 735 localPoints[nL] = localFace; 736 remotePoints[nL].rank = neighbor; 737 remotePoints[nL].index = remoteFace; 738 } 739 } 740 } else if (yp > 0) { /* top */ 741 if (zp < 0) { /* back */ 742 nleavesCheck += nVx; /* top back vertices */ 743 for (xv = 0; xv < nVx; ++xv, ++nL) { 744 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; 745 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 746 747 localPoints[nL] = localVertex; 748 remotePoints[nL].rank = neighbor; 749 remotePoints[nL].index = remoteVertex; 750 } 751 } else if (zp > 0) { /* front */ 752 nleavesCheck += nVx; /* top front vertices */ 753 for (xv = 0; xv < nVx; ++xv, ++nL) { 754 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; 755 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 756 757 localPoints[nL] = localVertex; 758 remotePoints[nL].rank = neighbor; 759 remotePoints[nL].index = remoteVertex; 760 } 761 } else { 762 nleavesCheck += nVx*nVz; /* top vertices */ 763 for (zv = 0; zv < nVz; ++zv) { 764 for (xv = 0; xv < nVx; ++xv, ++nL) { 765 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + xv + nC; 766 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 767 768 localPoints[nL] = localVertex; 769 remotePoints[nL].rank = neighbor; 770 remotePoints[nL].index = remoteVertex; 771 } 772 } 773 nleavesCheck += nyF; /* top faces */ 774 for (yf = 0; yf < nyF; ++yf, ++nL) { 775 /* THIS IS WRONG */ 776 const PetscInt localFace = 0 + nC+nV; 777 const PetscInt remoteFace = 0 + nC+nV; 778 779 localPoints[nL] = localFace; 780 remotePoints[nL].rank = neighbor; 781 remotePoints[nL].index = remoteFace; 782 } 783 } 784 } else { 785 if (zp < 0) { /* back */ 786 nleavesCheck += nVx*nVy; /* back vertices */ 787 for (yv = 0; yv < nVy; ++yv) { 788 for (xv = 0; xv < nVx; ++xv, ++nL) { 789 const PetscInt localVertex = ( 0*nVy + yv)*nVx + xv + nC; 790 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 791 792 localPoints[nL] = localVertex; 793 remotePoints[nL].rank = neighbor; 794 remotePoints[nL].index = remoteVertex; 795 } 796 } 797 nleavesCheck += nzF; /* back faces */ 798 for (zf = 0; zf < nzF; ++zf, ++nL) { 799 /* THIS IS WRONG */ 800 const PetscInt localFace = 0 + nC+nV; 801 const PetscInt remoteFace = 0 + nC+nV; 802 803 localPoints[nL] = localFace; 804 remotePoints[nL].rank = neighbor; 805 remotePoints[nL].index = remoteFace; 806 } 807 } else if (zp > 0) { /* front */ 808 nleavesCheck += nVx*nVy; /* front vertices */ 809 for (yv = 0; yv < nVy; ++yv) { 810 for (xv = 0; xv < nVx; ++xv, ++nL) { 811 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; 812 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 813 814 localPoints[nL] = localVertex; 815 remotePoints[nL].rank = neighbor; 816 remotePoints[nL].index = remoteVertex; 817 } 818 } 819 nleavesCheck += nzF; /* front faces */ 820 for (zf = 0; zf < nzF; ++zf, ++nL) { 821 /* THIS IS WRONG */ 822 const PetscInt localFace = 0 + nC+nV; 823 const PetscInt remoteFace = 0 + nC+nV; 824 825 localPoints[nL] = localFace; 826 remotePoints[nL].rank = neighbor; 827 remotePoints[nL].index = remoteFace; 828 } 829 } else { 830 /* Nothing is shared from the interior */ 831 } 832 } 833 } 834 } 835 } 836 } 837 } 838 /* TODO: Remove duplication in leaf determination */ 839 if (nleaves != nleavesCheck) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "The number of leaves %d did not match the number of remote leaves %d", nleaves, nleavesCheck); 840 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf);CHKERRQ(ierr); 841 ierr = PetscSFSetGraph(sf, pEnd, nleaves, localPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 842 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 843 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 844 *s = section; 845 PetscFunctionReturn(0); 846 } 847 848 #undef __FUNCT__ 849 #define __FUNCT__ "DMDASetVertexCoordinates" 850 PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu) 851 { 852 DM_DA *da = (DM_DA *) dm->data; 853 Vec coordinates; 854 PetscSection section; 855 PetscScalar *coords; 856 PetscReal h[3]; 857 PetscInt dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k; 858 PetscErrorCode ierr; 859 860 PetscFunctionBegin; 861 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 862 ierr = DMDAGetInfo(dm, &dim, &M, &N, &P, 0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 863 h[0] = (xu - xl)/M; 864 h[1] = (yu - yl)/N; 865 h[2] = (zu - zl)/P; 866 ierr = DMDAGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 867 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr); 868 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ion);CHKERRQ(ierr); 869 ierr = PetscSectionSetNumFields(section, 1);CHKERRQ(ierr); 870 ierr = PetscSectionSetFieldComponents(section, 0, dim);CHKERRQ(ierr); 871 ierr = PetscSectionSetChart(section, vStart, vEnd);CHKERRQ(ierr); 872 for (v = vStart; v < vEnd; ++v) { 873 ierr = PetscSectionSetDof(section, v, dim);CHKERRQ(ierr); 874 } 875 ierr = PetscSectionSetUp(section);CHKERRQ(ierr); 876 ierr = PetscSectionGetStorageSize(section, &size);CHKERRQ(ierr); 877 ierr = VecCreateSeq(PETSC_COMM_SELF, size, &coordinates);CHKERRQ(ierr); 878 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 879 for (k = 0; k < nVz; ++k) { 880 PetscInt ind[3], d, off; 881 882 ind[0] = 0; 883 ind[1] = 0; 884 ind[2] = k + da->zs; 885 for (j = 0; j < nVy; ++j) { 886 ind[1] = j + da->ys; 887 for (i = 0; i < nVx; ++i) { 888 const PetscInt vertex = (k*nVy + j)*nVx + i + vStart; 889 890 ierr = PetscSectionGetOffset(section, vertex, &off);CHKERRQ(ierr); 891 ind[0] = i + da->xs; 892 for (d = 0; d < dim; ++d) { 893 coords[off+d] = h[d]*ind[d]; 894 } 895 } 896 } 897 } 898 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 899 ierr = DMSetCoordinateSection(dm, section);CHKERRQ(ierr); 900 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 901 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr); 902 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 903 PetscFunctionReturn(0); 904 } 905 906 #undef __FUNCT__ 907 #define __FUNCT__ "DMDAProjectFunctionLocal" 908 PetscErrorCode DMDAProjectFunctionLocal(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 909 { 910 PetscDualSpace *sp; 911 PetscQuadrature q; 912 PetscSection section; 913 PetscScalar *values; 914 PetscReal *v0, *J, *detJ; 915 PetscInt numFields, numComp, dim, spDim, totDim = 0, numValues, cStart, cEnd, f, c, v, d; 916 PetscErrorCode ierr; 917 918 PetscFunctionBegin; 919 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 920 ierr = PetscFEGetQuadrature(fe[0], &q);CHKERRQ(ierr); 921 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 922 ierr = PetscMalloc(numFields * sizeof(PetscDualSpace), &sp);CHKERRQ(ierr); 923 for (f = 0; f < numFields; ++f) { 924 ierr = PetscFEGetDualSpace(fe[f], &sp[f]);CHKERRQ(ierr); 925 ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr); 926 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 927 totDim += spDim*numComp; 928 } 929 ierr = DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 930 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 931 ierr = DMDAVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr); 932 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim); 933 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 934 ierr = PetscMalloc3(dim*q.numPoints,&v0,dim*dim*q.numPoints,&J,q.numPoints,&detJ);CHKERRQ(ierr); 935 for (c = cStart; c < cEnd; ++c) { 936 PetscCellGeometry geom; 937 938 ierr = DMDAComputeCellGeometry(dm, c, &q, v0, J, NULL, detJ);CHKERRQ(ierr); 939 geom.v0 = v0; 940 geom.J = J; 941 geom.detJ = detJ; 942 for (f = 0, v = 0; f < numFields; ++f) { 943 void * const ctx = ctxs ? ctxs[f] : NULL; 944 ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr); 945 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 946 for (d = 0; d < spDim; ++d) { 947 ierr = PetscDualSpaceApply(sp[f], d, geom, numComp, funcs[f], ctx, &values[v]);CHKERRQ(ierr); 948 v += numComp; 949 } 950 } 951 ierr = DMDAVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr); 952 } 953 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 954 ierr = PetscFree3(v0,J,detJ);CHKERRQ(ierr); 955 ierr = PetscFree(sp);CHKERRQ(ierr); 956 PetscFunctionReturn(0); 957 } 958 959 #undef __FUNCT__ 960 #define __FUNCT__ "DMDAProjectFunction" 961 /*@C 962 DMDAProjectFunction - This projects the given function into the function space provided. 963 964 Input Parameters: 965 + dm - The DM 966 . fe - The PetscFE associated with the field 967 . funcs - The coordinate functions to evaluate 968 . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null. 969 - mode - The insertion mode for values 970 971 Output Parameter: 972 . X - vector 973 974 Level: developer 975 976 .seealso: DMDAComputeL2Diff() 977 @*/ 978 PetscErrorCode DMDAProjectFunction(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X) 979 { 980 Vec localX; 981 PetscErrorCode ierr; 982 983 PetscFunctionBegin; 984 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 985 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 986 ierr = DMDAProjectFunctionLocal(dm, fe, funcs, ctxs, mode, localX);CHKERRQ(ierr); 987 ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr); 988 ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr); 989 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 990 PetscFunctionReturn(0); 991 } 992 993 #undef __FUNCT__ 994 #define __FUNCT__ "DMDAComputeL2Diff" 995 /*@C 996 DMDAComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h. 997 998 Input Parameters: 999 + dm - The DM 1000 . fe - The PetscFE object for each field 1001 . funcs - The functions to evaluate for each field component 1002 . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null. 1003 - X - The coefficient vector u_h 1004 1005 Output Parameter: 1006 . diff - The diff ||u - u_h||_2 1007 1008 Level: developer 1009 1010 .seealso: DMDAProjectFunction() 1011 @*/ 1012 PetscErrorCode DMDAComputeL2Diff(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 1013 { 1014 const PetscInt debug = 0; 1015 PetscSection section; 1016 PetscQuadrature quad; 1017 Vec localX; 1018 PetscScalar *funcVal; 1019 PetscReal *coords, *v0, *J, *invJ, detJ; 1020 PetscReal localDiff = 0.0; 1021 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp; 1022 PetscErrorCode ierr; 1023 1024 PetscFunctionBegin; 1025 ierr = DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 1026 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1027 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1028 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1029 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1030 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1031 for (field = 0; field < numFields; ++field) { 1032 PetscInt Nc; 1033 1034 ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr); 1035 numComponents += Nc; 1036 } 1037 /* There are no BC values in DAs right now: ierr = DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 1038 ierr = PetscMalloc5(numComponents,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1039 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1040 ierr = PetscFEGetQuadrature(fe[0], &quad);CHKERRQ(ierr); 1041 for (c = cStart; c < cEnd; ++c) { 1042 PetscScalar *x = NULL; 1043 PetscReal elemDiff = 0.0; 1044 1045 ierr = DMDAComputeCellGeometry(dm, c, &quad, v0, J, invJ, &detJ);CHKERRQ(ierr); 1046 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 1047 ierr = DMDAVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1048 1049 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 1050 void * const ctx = ctxs ? ctxs[field] : NULL; 1051 const PetscInt numQuadPoints = quad.numPoints; 1052 const PetscReal *quadPoints = quad.points; 1053 const PetscReal *quadWeights = quad.weights; 1054 PetscReal *basis; 1055 PetscInt numBasisFuncs, numBasisComps, q, d, e, fc, f; 1056 1057 ierr = PetscFEGetDimension(fe[field], &numBasisFuncs);CHKERRQ(ierr); 1058 ierr = PetscFEGetNumComponents(fe[field], &numBasisComps);CHKERRQ(ierr); 1059 ierr = PetscFEGetDefaultTabulation(fe[field], &basis, NULL, NULL);CHKERRQ(ierr); 1060 if (debug) { 1061 char title[1024]; 1062 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 1063 ierr = DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);CHKERRQ(ierr); 1064 } 1065 for (q = 0; q < numQuadPoints; ++q) { 1066 for (d = 0; d < dim; d++) { 1067 coords[d] = v0[d]; 1068 for (e = 0; e < dim; e++) { 1069 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 1070 } 1071 } 1072 (*funcs[field])(coords, funcVal, ctx); 1073 for (fc = 0; fc < numBasisComps; ++fc) { 1074 PetscScalar interpolant = 0.0; 1075 1076 for (f = 0; f < numBasisFuncs; ++f) { 1077 const PetscInt fidx = f*numBasisComps+fc; 1078 interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx]; 1079 } 1080 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);} 1081 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 1082 } 1083 } 1084 comp += numBasisComps; 1085 fieldOffset += numBasisFuncs*numBasisComps; 1086 } 1087 ierr = DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1088 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 1089 localDiff += elemDiff; 1090 } 1091 ierr = PetscFree5(funcVal,coords,v0,J,invJ);CHKERRQ(ierr); 1092 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1093 ierr = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1094 *diff = PetscSqrtReal(*diff); 1095 PetscFunctionReturn(0); 1096 } 1097 1098 /* ------------------------------------------------------------------- */ 1099 1100 #undef __FUNCT__ 1101 #define __FUNCT__ "DMDAGetArray" 1102 /*@C 1103 DMDAGetArray - Gets a work array for a DMDA 1104 1105 Input Parameter: 1106 + da - information about my local patch 1107 - ghosted - do you want arrays for the ghosted or nonghosted patch 1108 1109 Output Parameters: 1110 . vptr - array data structured 1111 1112 Note: The vector values are NOT initialized and may have garbage in them, so you may need 1113 to zero them. 1114 1115 Level: advanced 1116 1117 .seealso: DMDARestoreArray() 1118 1119 @*/ 1120 PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr) 1121 { 1122 PetscErrorCode ierr; 1123 PetscInt j,i,xs,ys,xm,ym,zs,zm; 1124 char *iarray_start; 1125 void **iptr = (void**)vptr; 1126 DM_DA *dd = (DM_DA*)da->data; 1127 1128 PetscFunctionBegin; 1129 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1130 if (ghosted) { 1131 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1132 if (dd->arrayghostedin[i]) { 1133 *iptr = dd->arrayghostedin[i]; 1134 iarray_start = (char*)dd->startghostedin[i]; 1135 dd->arrayghostedin[i] = NULL; 1136 dd->startghostedin[i] = NULL; 1137 1138 goto done; 1139 } 1140 } 1141 xs = dd->Xs; 1142 ys = dd->Ys; 1143 zs = dd->Zs; 1144 xm = dd->Xe-dd->Xs; 1145 ym = dd->Ye-dd->Ys; 1146 zm = dd->Ze-dd->Zs; 1147 } else { 1148 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1149 if (dd->arrayin[i]) { 1150 *iptr = dd->arrayin[i]; 1151 iarray_start = (char*)dd->startin[i]; 1152 dd->arrayin[i] = NULL; 1153 dd->startin[i] = NULL; 1154 1155 goto done; 1156 } 1157 } 1158 xs = dd->xs; 1159 ys = dd->ys; 1160 zs = dd->zs; 1161 xm = dd->xe-dd->xs; 1162 ym = dd->ye-dd->ys; 1163 zm = dd->ze-dd->zs; 1164 } 1165 1166 switch (dd->dim) { 1167 case 1: { 1168 void *ptr; 1169 1170 ierr = PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1171 1172 ptr = (void*)(iarray_start - xs*sizeof(PetscScalar)); 1173 *iptr = (void*)ptr; 1174 break; 1175 } 1176 case 2: { 1177 void **ptr; 1178 1179 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1180 1181 ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*)); 1182 for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs); 1183 *iptr = (void*)ptr; 1184 break; 1185 } 1186 case 3: { 1187 void ***ptr,**bptr; 1188 1189 ierr = PetscMalloc((zm+1)*sizeof(void**)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1190 1191 ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*)); 1192 bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**)); 1193 for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys); 1194 for (i=zs; i<zs+zm; i++) { 1195 for (j=ys; j<ys+ym; j++) { 1196 ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs); 1197 } 1198 } 1199 *iptr = (void*)ptr; 1200 break; 1201 } 1202 default: 1203 SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 1204 } 1205 1206 done: 1207 /* add arrays to the checked out list */ 1208 if (ghosted) { 1209 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1210 if (!dd->arrayghostedout[i]) { 1211 dd->arrayghostedout[i] = *iptr; 1212 dd->startghostedout[i] = iarray_start; 1213 break; 1214 } 1215 } 1216 } else { 1217 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1218 if (!dd->arrayout[i]) { 1219 dd->arrayout[i] = *iptr; 1220 dd->startout[i] = iarray_start; 1221 break; 1222 } 1223 } 1224 } 1225 PetscFunctionReturn(0); 1226 } 1227 1228 #undef __FUNCT__ 1229 #define __FUNCT__ "DMDARestoreArray" 1230 /*@C 1231 DMDARestoreArray - Restores an array of derivative types for a DMDA 1232 1233 Input Parameter: 1234 + da - information about my local patch 1235 . ghosted - do you want arrays for the ghosted or nonghosted patch 1236 - vptr - array data structured to be passed to ad_FormFunctionLocal() 1237 1238 Level: advanced 1239 1240 .seealso: DMDAGetArray() 1241 1242 @*/ 1243 PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr) 1244 { 1245 PetscInt i; 1246 void **iptr = (void**)vptr,*iarray_start = 0; 1247 DM_DA *dd = (DM_DA*)da->data; 1248 1249 PetscFunctionBegin; 1250 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1251 if (ghosted) { 1252 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1253 if (dd->arrayghostedout[i] == *iptr) { 1254 iarray_start = dd->startghostedout[i]; 1255 dd->arrayghostedout[i] = NULL; 1256 dd->startghostedout[i] = NULL; 1257 break; 1258 } 1259 } 1260 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1261 if (!dd->arrayghostedin[i]) { 1262 dd->arrayghostedin[i] = *iptr; 1263 dd->startghostedin[i] = iarray_start; 1264 break; 1265 } 1266 } 1267 } else { 1268 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1269 if (dd->arrayout[i] == *iptr) { 1270 iarray_start = dd->startout[i]; 1271 dd->arrayout[i] = NULL; 1272 dd->startout[i] = NULL; 1273 break; 1274 } 1275 } 1276 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1277 if (!dd->arrayin[i]) { 1278 dd->arrayin[i] = *iptr; 1279 dd->startin[i] = iarray_start; 1280 break; 1281 } 1282 } 1283 } 1284 PetscFunctionReturn(0); 1285 } 1286 1287