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