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 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 PetscErrorCode DMCreateLocalVector_DA(DM da,Vec *g) 49 { 50 PetscErrorCode ierr; 51 DM_DA *dd = (DM_DA*)da->data; 52 53 PetscFunctionBegin; 54 PetscValidHeaderSpecific(da,DM_CLASSID,1); 55 PetscValidPointer(g,2); 56 if (da->defaultSection) { 57 ierr = DMCreateLocalVector_Section_Private(da,g);CHKERRQ(ierr); 58 } else { 59 ierr = VecCreate(PETSC_COMM_SELF,g);CHKERRQ(ierr); 60 ierr = VecSetSizes(*g,dd->nlocal,PETSC_DETERMINE);CHKERRQ(ierr); 61 ierr = VecSetBlockSize(*g,dd->w);CHKERRQ(ierr); 62 ierr = VecSetType(*g,da->vectype);CHKERRQ(ierr); 63 ierr = VecSetDM(*g, da);CHKERRQ(ierr); 64 #if defined(PETSC_HAVE_MATLAB_ENGINE) 65 if (dd->w == 1 && da->dim == 2) { 66 ierr = PetscObjectComposeFunction((PetscObject)*g,"PetscMatlabEnginePut_C",VecMatlabEnginePut_DA2d);CHKERRQ(ierr); 67 } 68 #endif 69 } 70 PetscFunctionReturn(0); 71 } 72 73 /*@ 74 DMDAGetNumCells - Get the number of cells in the local piece of the DMDA. This includes ghost cells. 75 76 Input Parameter: 77 . dm - The DM object 78 79 Output Parameters: 80 + numCellsX - The number of local cells in the x-direction 81 . numCellsY - The number of local cells in the y-direction 82 . numCellsZ - The number of local cells in the z-direction 83 - numCells - The number of local cells 84 85 Level: developer 86 87 .seealso: DMDAGetCellPoint() 88 @*/ 89 PetscErrorCode DMDAGetNumCells(DM dm, PetscInt *numCellsX, PetscInt *numCellsY, PetscInt *numCellsZ, PetscInt *numCells) 90 { 91 DM_DA *da = (DM_DA*) dm->data; 92 const PetscInt dim = dm->dim; 93 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 94 const PetscInt nC = (mx)*(dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1); 95 96 PetscFunctionBegin; 97 PetscValidHeaderSpecificType(dm, DM_CLASSID, 1,DMDA); 98 if (numCellsX) { 99 PetscValidIntPointer(numCellsX,2); 100 *numCellsX = mx; 101 } 102 if (numCellsY) { 103 PetscValidIntPointer(numCellsY,3); 104 *numCellsY = my; 105 } 106 if (numCellsZ) { 107 PetscValidIntPointer(numCellsZ,4); 108 *numCellsZ = mz; 109 } 110 if (numCells) { 111 PetscValidIntPointer(numCells,5); 112 *numCells = nC; 113 } 114 PetscFunctionReturn(0); 115 } 116 117 /*@ 118 DMDAGetCellPoint - Get the DM point corresponding to the tuple (i, j, k) in the DMDA 119 120 Input Parameters: 121 + dm - The DM object 122 - i,j,k - The global indices for the cell 123 124 Output Parameters: 125 . point - The local DM point 126 127 Level: developer 128 129 .seealso: DMDAGetNumCells() 130 @*/ 131 PetscErrorCode DMDAGetCellPoint(DM dm, PetscInt i, PetscInt j, PetscInt k, PetscInt *point) 132 { 133 const PetscInt dim = dm->dim; 134 DMDALocalInfo info; 135 PetscErrorCode ierr; 136 137 PetscFunctionBegin; 138 PetscValidHeaderSpecificType(dm, DM_CLASSID, 1,DMDA); 139 PetscValidIntPointer(point,5); 140 ierr = DMDAGetLocalInfo(dm, &info);CHKERRQ(ierr); 141 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);} 142 if (dim > 1) {if ((j < info.gys) || (j >= info.gys+info.gym)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Y index %d not in [%d, %d)", j, info.gys, info.gys+info.gym);} 143 if (dim > 2) {if ((k < info.gzs) || (k >= info.gzs+info.gzm)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Z index %d not in [%d, %d)", k, info.gzs, info.gzs+info.gzm);} 144 *point = i + (dim > 1 ? (j + (dim > 2 ? k*info.gym : 0))*info.gxm : 0); 145 PetscFunctionReturn(0); 146 } 147 148 PetscErrorCode DMDAGetNumVertices(DM dm, PetscInt *numVerticesX, PetscInt *numVerticesY, PetscInt *numVerticesZ, PetscInt *numVertices) 149 { 150 DM_DA *da = (DM_DA*) dm->data; 151 const PetscInt dim = dm->dim; 152 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 153 const PetscInt nVx = mx+1; 154 const PetscInt nVy = dim > 1 ? (my+1) : 1; 155 const PetscInt nVz = dim > 2 ? (mz+1) : 1; 156 const PetscInt nV = nVx*nVy*nVz; 157 158 PetscFunctionBegin; 159 if (numVerticesX) { 160 PetscValidIntPointer(numVerticesX,2); 161 *numVerticesX = nVx; 162 } 163 if (numVerticesY) { 164 PetscValidIntPointer(numVerticesY,3); 165 *numVerticesY = nVy; 166 } 167 if (numVerticesZ) { 168 PetscValidIntPointer(numVerticesZ,4); 169 *numVerticesZ = nVz; 170 } 171 if (numVertices) { 172 PetscValidIntPointer(numVertices,5); 173 *numVertices = nV; 174 } 175 PetscFunctionReturn(0); 176 } 177 178 PetscErrorCode DMDAGetNumFaces(DM dm, PetscInt *numXFacesX, PetscInt *numXFaces, PetscInt *numYFacesY, PetscInt *numYFaces, PetscInt *numZFacesZ, PetscInt *numZFaces) 179 { 180 DM_DA *da = (DM_DA*) dm->data; 181 const PetscInt dim = dm->dim; 182 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 183 const PetscInt nxF = (dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1); 184 const PetscInt nXF = (mx+1)*nxF; 185 const PetscInt nyF = mx*(dim > 2 ? mz : 1); 186 const PetscInt nYF = dim > 1 ? (my+1)*nyF : 0; 187 const PetscInt nzF = mx*(dim > 1 ? my : 0); 188 const PetscInt nZF = dim > 2 ? (mz+1)*nzF : 0; 189 190 PetscFunctionBegin; 191 if (numXFacesX) { 192 PetscValidIntPointer(numXFacesX,2); 193 *numXFacesX = nxF; 194 } 195 if (numXFaces) { 196 PetscValidIntPointer(numXFaces,3); 197 *numXFaces = nXF; 198 } 199 if (numYFacesY) { 200 PetscValidIntPointer(numYFacesY,4); 201 *numYFacesY = nyF; 202 } 203 if (numYFaces) { 204 PetscValidIntPointer(numYFaces,5); 205 *numYFaces = nYF; 206 } 207 if (numZFacesZ) { 208 PetscValidIntPointer(numZFacesZ,6); 209 *numZFacesZ = nzF; 210 } 211 if (numZFaces) { 212 PetscValidIntPointer(numZFaces,7); 213 *numZFaces = nZF; 214 } 215 PetscFunctionReturn(0); 216 } 217 218 PetscErrorCode DMDAGetHeightStratum(DM dm, PetscInt height, PetscInt *pStart, PetscInt *pEnd) 219 { 220 const PetscInt dim = dm->dim; 221 PetscInt nC, nV, nXF, nYF, nZF; 222 PetscErrorCode ierr; 223 224 PetscFunctionBegin; 225 if (pStart) PetscValidIntPointer(pStart,3); 226 if (pEnd) PetscValidIntPointer(pEnd,4); 227 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); 228 ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr); 229 ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr); 230 if (height == 0) { 231 /* Cells */ 232 if (pStart) *pStart = 0; 233 if (pEnd) *pEnd = nC; 234 } else if (height == 1) { 235 /* Faces */ 236 if (pStart) *pStart = nC+nV; 237 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 238 } else if (height == dim) { 239 /* Vertices */ 240 if (pStart) *pStart = nC; 241 if (pEnd) *pEnd = nC+nV; 242 } else if (height < 0) { 243 /* All points */ 244 if (pStart) *pStart = 0; 245 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 246 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of height %d in the DA", height); 247 PetscFunctionReturn(0); 248 } 249 250 PetscErrorCode DMDAGetDepthStratum(DM dm, PetscInt depth, PetscInt *pStart, PetscInt *pEnd) 251 { 252 const PetscInt dim = dm->dim; 253 PetscInt nC, nV, nXF, nYF, nZF; 254 PetscErrorCode ierr; 255 256 PetscFunctionBegin; 257 if (pStart) PetscValidIntPointer(pStart,3); 258 if (pEnd) PetscValidIntPointer(pEnd,4); 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 if (depth == dim) { 263 /* Cells */ 264 if (pStart) *pStart = 0; 265 if (pEnd) *pEnd = nC; 266 } else if (depth == dim-1) { 267 /* Faces */ 268 if (pStart) *pStart = nC+nV; 269 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 270 } else if (depth == 0) { 271 /* Vertices */ 272 if (pStart) *pStart = nC; 273 if (pEnd) *pEnd = nC+nV; 274 } else if (depth < 0) { 275 /* All points */ 276 if (pStart) *pStart = 0; 277 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 278 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of depth %d in the DA", depth); 279 PetscFunctionReturn(0); 280 } 281 282 PetscErrorCode DMDAGetConeSize(DM dm, PetscInt p, PetscInt *coneSize) 283 { 284 const PetscInt dim = dm->dim; 285 PetscInt nC, nV, nXF, nYF, nZF; 286 PetscErrorCode ierr; 287 288 PetscFunctionBegin; 289 *coneSize = 0; 290 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); 291 ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr); 292 ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr); 293 switch (dim) { 294 case 2: 295 if (p >= 0) { 296 if (p < nC) { 297 *coneSize = 4; 298 } else if (p < nC+nV) { 299 *coneSize = 0; 300 } else if (p < nC+nV+nXF+nYF+nZF) { 301 *coneSize = 2; 302 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF); 303 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p); 304 break; 305 case 3: 306 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D"); 307 break; 308 } 309 PetscFunctionReturn(0); 310 } 311 312 PetscErrorCode DMDAGetCone(DM dm, PetscInt p, PetscInt *cone[]) 313 { 314 const PetscInt dim = dm->dim; 315 PetscInt nCx, nCy, nCz, nC, nVx, nVy, nVz, nV, nxF, nyF, nzF, nXF, nYF, nZF; 316 PetscErrorCode ierr; 317 318 PetscFunctionBegin; 319 if (!cone) {ierr = DMGetWorkArray(dm, 6, MPIU_INT, cone);CHKERRQ(ierr);} 320 ierr = DMDAGetNumCells(dm, &nCx, &nCy, &nCz, &nC);CHKERRQ(ierr); 321 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr); 322 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 323 switch (dim) { 324 case 2: 325 if (p >= 0) { 326 if (p < nC) { 327 const PetscInt cy = p / nCx; 328 const PetscInt cx = p % nCx; 329 330 (*cone)[0] = cy*nxF + cx + nC+nV; 331 (*cone)[1] = cx*nyF + cy + nyF + nC+nV+nXF; 332 (*cone)[2] = cy*nxF + cx + nxF + nC+nV; 333 (*cone)[3] = cx*nyF + cy + nC+nV+nXF; 334 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do cell cones"); 335 } else if (p < nC+nV) { 336 } else if (p < nC+nV+nXF) { 337 const PetscInt fy = (p - nC+nV) / nxF; 338 const PetscInt fx = (p - nC+nV) % nxF; 339 340 (*cone)[0] = fy*nVx + fx + nC; 341 (*cone)[1] = fy*nVx + fx + 1 + nC; 342 } else if (p < nC+nV+nXF+nYF) { 343 const PetscInt fx = (p - nC+nV+nXF) / nyF; 344 const PetscInt fy = (p - nC+nV+nXF) % nyF; 345 346 (*cone)[0] = fy*nVx + fx + nC; 347 (*cone)[1] = fy*nVx + fx + nVx + nC; 348 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF); 349 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p); 350 break; 351 case 3: 352 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D"); 353 break; 354 } 355 PetscFunctionReturn(0); 356 } 357 358 PetscErrorCode DMDARestoreCone(DM dm, PetscInt p, PetscInt *cone[]) 359 { 360 PetscErrorCode ierr; 361 362 PetscFunctionBegin; 363 ierr = DMGetWorkArray(dm, 6, MPIU_INT, cone);CHKERRQ(ierr); 364 PetscFunctionReturn(0); 365 } 366 367 PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu) 368 { 369 DM_DA *da = (DM_DA *) dm->data; 370 Vec coordinates; 371 PetscSection section; 372 PetscScalar *coords; 373 PetscReal h[3]; 374 PetscInt dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k; 375 PetscErrorCode ierr; 376 377 PetscFunctionBegin; 378 PetscValidHeaderSpecificType(dm, DM_CLASSID, 1,DMDA); 379 ierr = DMDAGetInfo(dm, &dim, &M, &N, &P, 0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 380 if (dim > 3) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"The following code only works for dim <= 3"); 381 h[0] = (xu - xl)/M; 382 h[1] = (yu - yl)/N; 383 h[2] = (zu - zl)/P; 384 ierr = DMDAGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 385 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr); 386 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ion);CHKERRQ(ierr); 387 ierr = PetscSectionSetNumFields(section, 1);CHKERRQ(ierr); 388 ierr = PetscSectionSetFieldComponents(section, 0, dim);CHKERRQ(ierr); 389 ierr = PetscSectionSetChart(section, vStart, vEnd);CHKERRQ(ierr); 390 for (v = vStart; v < vEnd; ++v) { 391 ierr = PetscSectionSetDof(section, v, dim);CHKERRQ(ierr); 392 } 393 ierr = PetscSectionSetUp(section);CHKERRQ(ierr); 394 ierr = PetscSectionGetStorageSize(section, &size);CHKERRQ(ierr); 395 ierr = VecCreateSeq(PETSC_COMM_SELF, size, &coordinates);CHKERRQ(ierr); 396 ierr = PetscObjectSetName((PetscObject)coordinates,"coordinates");CHKERRQ(ierr); 397 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 398 for (k = 0; k < nVz; ++k) { 399 PetscInt ind[3], d, off; 400 401 ind[0] = 0; 402 ind[1] = 0; 403 ind[2] = k + da->zs; 404 for (j = 0; j < nVy; ++j) { 405 ind[1] = j + da->ys; 406 for (i = 0; i < nVx; ++i) { 407 const PetscInt vertex = (k*nVy + j)*nVx + i + vStart; 408 409 ierr = PetscSectionGetOffset(section, vertex, &off);CHKERRQ(ierr); 410 ind[0] = i + da->xs; 411 for (d = 0; d < dim; ++d) { 412 coords[off+d] = h[d]*ind[d]; 413 } 414 } 415 } 416 } 417 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 418 ierr = DMSetCoordinateSection(dm, PETSC_DETERMINE, section);CHKERRQ(ierr); 419 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 420 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr); 421 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 422 PetscFunctionReturn(0); 423 } 424 425 PetscErrorCode DMProjectFunctionLocal_DA(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 426 { 427 PetscDS prob; 428 PetscFE fe; 429 PetscDualSpace sp; 430 PetscQuadrature q; 431 PetscSection section; 432 PetscScalar *values; 433 PetscInt numFields, Nc, dim, dimEmbed, spDim, totDim, numValues, cStart, cEnd, f, c, v, d; 434 PetscFEGeom *geom; 435 PetscErrorCode ierr; 436 437 PetscFunctionBegin; 438 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 439 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 440 ierr = PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);CHKERRQ(ierr); 441 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 442 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 443 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 444 ierr = DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 445 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 446 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 447 ierr = DMDAVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr); 448 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim); 449 ierr = DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr); 450 ierr = PetscFEGeomCreate(q,1,dimEmbed,PETSC_FALSE,&geom);CHKERRQ(ierr); 451 for (c = cStart; c < cEnd; ++c) { 452 ierr = DMDAComputeCellGeometryFEM(dm, c, q, geom->v, geom->J, NULL, geom->detJ);CHKERRQ(ierr); 453 for (f = 0, v = 0; f < numFields; ++f) { 454 void * const ctx = ctxs ? ctxs[f] : NULL; 455 456 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 457 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 458 ierr = PetscFEGetDualSpace(fe, &sp);CHKERRQ(ierr); 459 ierr = PetscDualSpaceGetDimension(sp, &spDim);CHKERRQ(ierr); 460 for (d = 0; d < spDim; ++d) { 461 ierr = PetscDualSpaceApply(sp, d, time, geom, Nc, funcs[f], ctx, &values[v]);CHKERRQ(ierr); 462 } 463 } 464 ierr = DMDAVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr); 465 } 466 ierr = PetscFEGeomDestroy(&geom);CHKERRQ(ierr); 467 ierr = DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);CHKERRQ(ierr); 468 PetscFunctionReturn(0); 469 } 470 471 PetscErrorCode DMComputeL2Diff_DA(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 472 { 473 const PetscInt debug = 0; 474 PetscDS prob; 475 PetscFE fe; 476 PetscQuadrature quad; 477 PetscSection section; 478 Vec localX; 479 PetscScalar *funcVal; 480 PetscReal *coords, *v0, *J, *invJ, detJ; 481 PetscReal localDiff = 0.0; 482 PetscInt dim, numFields, Nc, cStart, cEnd, c, field, fieldOffset, comp; 483 PetscErrorCode ierr; 484 485 PetscFunctionBegin; 486 ierr = DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 487 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 488 ierr = PetscDSGetTotalComponents(prob, &Nc);CHKERRQ(ierr); 489 ierr = PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);CHKERRQ(ierr); 490 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 491 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 492 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 493 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 494 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 495 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 496 /* There are no BC values in DAs right now: ierr = DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 497 ierr = PetscMalloc5(Nc,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 498 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 499 for (c = cStart; c < cEnd; ++c) { 500 PetscScalar *x = NULL; 501 PetscReal elemDiff = 0.0; 502 503 ierr = DMDAComputeCellGeometryFEM(dm, c, quad, v0, J, invJ, &detJ);CHKERRQ(ierr); 504 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 505 ierr = DMDAVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 506 507 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 508 void * const ctx = ctxs ? ctxs[field] : NULL; 509 const PetscReal *quadPoints, *quadWeights; 510 PetscReal *basis; 511 PetscInt qNc, Nq, Nb, Nc, q, d, e, fc, f; 512 513 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 514 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 515 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 516 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 517 if (qNc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Quadrature components %D != %D field components\n", qNc, Nc); 518 ierr = PetscFEGetDefaultTabulation(fe, &basis, NULL, NULL);CHKERRQ(ierr); 519 if (debug) { 520 char title[1024]; 521 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 522 ierr = DMPrintCellVector(c, title, Nb, &x[fieldOffset]);CHKERRQ(ierr); 523 } 524 for (q = 0; q < Nq; ++q) { 525 for (d = 0; d < dim; d++) { 526 coords[d] = v0[d]; 527 for (e = 0; e < dim; e++) { 528 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 529 } 530 } 531 ierr = (*funcs[field])(dim, time, coords, numFields, funcVal, ctx);CHKERRQ(ierr); 532 for (fc = 0; fc < Nc; ++fc) { 533 PetscScalar interpolant = 0.0; 534 535 for (f = 0; f < Nb; ++f) { 536 interpolant += x[fieldOffset+f]*basis[(q*Nb+f)*Nc+fc]; 537 } 538 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q*Nc+fc]*detJ);CHKERRQ(ierr);} 539 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q*Nc+c]*detJ; 540 } 541 } 542 comp += Nc; 543 fieldOffset += Nb; 544 } 545 ierr = DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 546 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 547 localDiff += elemDiff; 548 } 549 ierr = PetscFree5(funcVal,coords,v0,J,invJ);CHKERRQ(ierr); 550 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 551 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 552 *diff = PetscSqrtReal(*diff); 553 PetscFunctionReturn(0); 554 } 555 556 PetscErrorCode DMComputeL2GradientDiff_DA(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff) 557 { 558 const PetscInt debug = 0; 559 PetscDS prob; 560 PetscFE fe; 561 PetscQuadrature quad; 562 PetscSection section; 563 Vec localX; 564 PetscScalar *funcVal, *interpolantVec; 565 PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ; 566 PetscReal localDiff = 0.0; 567 PetscInt dim, numFields, Nc, cStart, cEnd, c, field, fieldOffset, comp; 568 PetscErrorCode ierr; 569 570 PetscFunctionBegin; 571 ierr = DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 572 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 573 ierr = PetscDSGetTotalComponents(prob, &Nc);CHKERRQ(ierr); 574 ierr = PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);CHKERRQ(ierr); 575 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 576 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 577 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 578 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 579 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 580 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 581 /* There are no BC values in DAs right now: ierr = DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 582 ierr = PetscMalloc7(Nc,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr); 583 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 584 for (c = cStart; c < cEnd; ++c) { 585 PetscScalar *x = NULL; 586 PetscReal elemDiff = 0.0; 587 588 ierr = DMDAComputeCellGeometryFEM(dm, c, quad, v0, J, invJ, &detJ);CHKERRQ(ierr); 589 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 590 ierr = DMDAVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 591 592 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 593 void * const ctx = ctxs ? ctxs[field] : NULL; 594 const PetscReal *quadPoints, *quadWeights; 595 PetscReal *basisDer; 596 PetscInt qNc, Nq, Nb, Nc, q, d, e, fc, f, g; 597 598 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 599 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 600 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 601 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 602 if (qNc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Quadrature components %D != %D field components\n", qNc, Nc); 603 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 604 if (debug) { 605 char title[1024]; 606 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 607 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 608 } 609 for (q = 0; q < Nq; ++q) { 610 for (d = 0; d < dim; d++) { 611 coords[d] = v0[d]; 612 for (e = 0; e < dim; e++) { 613 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 614 } 615 } 616 ierr = (*funcs[field])(dim, time, coords, n, numFields, funcVal, ctx);CHKERRQ(ierr); 617 for (fc = 0; fc < Nc; ++fc) { 618 PetscScalar interpolant = 0.0; 619 620 for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0; 621 for (f = 0; f < Nb; ++f) { 622 for (d = 0; d < dim; ++d) { 623 realSpaceDer[d] = 0.0; 624 for (g = 0; g < dim; ++g) { 625 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb+f)*dim+g]; 626 } 627 interpolantVec[d] += x[fieldOffset+f]*realSpaceDer[d]; 628 } 629 } 630 for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d]; 631 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q*Nc+c]*detJ);CHKERRQ(ierr);} 632 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q*Nc+c]*detJ; 633 } 634 } 635 comp += Nc; 636 fieldOffset += Nb*Nc; 637 } 638 ierr = DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 639 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 640 localDiff += elemDiff; 641 } 642 ierr = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr); 643 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 644 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 645 *diff = PetscSqrtReal(*diff); 646 PetscFunctionReturn(0); 647 } 648 649 /* ------------------------------------------------------------------- */ 650 651 /*@C 652 DMDAGetArray - Gets a work array for a DMDA 653 654 Input Parameter: 655 + da - information about my local patch 656 - ghosted - do you want arrays for the ghosted or nonghosted patch 657 658 Output Parameters: 659 . vptr - array data structured 660 661 Note: The vector values are NOT initialized and may have garbage in them, so you may need 662 to zero them. 663 664 Level: advanced 665 666 .seealso: DMDARestoreArray() 667 668 @*/ 669 PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr) 670 { 671 PetscErrorCode ierr; 672 PetscInt j,i,xs,ys,xm,ym,zs,zm; 673 char *iarray_start; 674 void **iptr = (void**)vptr; 675 DM_DA *dd = (DM_DA*)da->data; 676 677 PetscFunctionBegin; 678 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA); 679 if (ghosted) { 680 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 681 if (dd->arrayghostedin[i]) { 682 *iptr = dd->arrayghostedin[i]; 683 iarray_start = (char*)dd->startghostedin[i]; 684 dd->arrayghostedin[i] = NULL; 685 dd->startghostedin[i] = NULL; 686 687 goto done; 688 } 689 } 690 xs = dd->Xs; 691 ys = dd->Ys; 692 zs = dd->Zs; 693 xm = dd->Xe-dd->Xs; 694 ym = dd->Ye-dd->Ys; 695 zm = dd->Ze-dd->Zs; 696 } else { 697 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 698 if (dd->arrayin[i]) { 699 *iptr = dd->arrayin[i]; 700 iarray_start = (char*)dd->startin[i]; 701 dd->arrayin[i] = NULL; 702 dd->startin[i] = NULL; 703 704 goto done; 705 } 706 } 707 xs = dd->xs; 708 ys = dd->ys; 709 zs = dd->zs; 710 xm = dd->xe-dd->xs; 711 ym = dd->ye-dd->ys; 712 zm = dd->ze-dd->zs; 713 } 714 715 switch (da->dim) { 716 case 1: { 717 void *ptr; 718 719 ierr = PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 720 721 ptr = (void*)(iarray_start - xs*sizeof(PetscScalar)); 722 *iptr = (void*)ptr; 723 break; 724 } 725 case 2: { 726 void **ptr; 727 728 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 729 730 ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*)); 731 for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs); 732 *iptr = (void*)ptr; 733 break; 734 } 735 case 3: { 736 void ***ptr,**bptr; 737 738 ierr = PetscMalloc((zm+1)*sizeof(void**)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 739 740 ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*)); 741 bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**)); 742 for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys); 743 for (i=zs; i<zs+zm; i++) { 744 for (j=ys; j<ys+ym; j++) { 745 ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs); 746 } 747 } 748 *iptr = (void*)ptr; 749 break; 750 } 751 default: 752 SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",da->dim); 753 } 754 755 done: 756 /* add arrays to the checked out list */ 757 if (ghosted) { 758 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 759 if (!dd->arrayghostedout[i]) { 760 dd->arrayghostedout[i] = *iptr; 761 dd->startghostedout[i] = iarray_start; 762 break; 763 } 764 } 765 } else { 766 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 767 if (!dd->arrayout[i]) { 768 dd->arrayout[i] = *iptr; 769 dd->startout[i] = iarray_start; 770 break; 771 } 772 } 773 } 774 PetscFunctionReturn(0); 775 } 776 777 /*@C 778 DMDARestoreArray - Restores an array of derivative types for a DMDA 779 780 Input Parameter: 781 + da - information about my local patch 782 . ghosted - do you want arrays for the ghosted or nonghosted patch 783 - vptr - array data structured to be passed to ad_FormFunctionLocal() 784 785 Level: advanced 786 787 .seealso: DMDAGetArray() 788 789 @*/ 790 PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr) 791 { 792 PetscInt i; 793 void **iptr = (void**)vptr,*iarray_start = 0; 794 DM_DA *dd = (DM_DA*)da->data; 795 796 PetscFunctionBegin; 797 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA); 798 if (ghosted) { 799 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 800 if (dd->arrayghostedout[i] == *iptr) { 801 iarray_start = dd->startghostedout[i]; 802 dd->arrayghostedout[i] = NULL; 803 dd->startghostedout[i] = NULL; 804 break; 805 } 806 } 807 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 808 if (!dd->arrayghostedin[i]) { 809 dd->arrayghostedin[i] = *iptr; 810 dd->startghostedin[i] = iarray_start; 811 break; 812 } 813 } 814 } else { 815 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 816 if (dd->arrayout[i] == *iptr) { 817 iarray_start = dd->startout[i]; 818 dd->arrayout[i] = NULL; 819 dd->startout[i] = NULL; 820 break; 821 } 822 } 823 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 824 if (!dd->arrayin[i]) { 825 dd->arrayin[i] = *iptr; 826 dd->startin[i] = iarray_start; 827 break; 828 } 829 } 830 } 831 PetscFunctionReturn(0); 832 } 833 834