1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petscsf.h> 3 4 #include <petsc/private/petscfeimpl.h> 5 #include <petsc/private/petscfvimpl.h> 6 7 #undef __FUNCT__ 8 #define __FUNCT__ "DMPlexGetScale" 9 PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale) 10 { 11 DM_Plex *mesh = (DM_Plex*) dm->data; 12 13 PetscFunctionBegin; 14 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 15 PetscValidPointer(scale, 3); 16 *scale = mesh->scale[unit]; 17 PetscFunctionReturn(0); 18 } 19 20 #undef __FUNCT__ 21 #define __FUNCT__ "DMPlexSetScale" 22 PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale) 23 { 24 DM_Plex *mesh = (DM_Plex*) dm->data; 25 26 PetscFunctionBegin; 27 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 28 mesh->scale[unit] = scale; 29 PetscFunctionReturn(0); 30 } 31 32 PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k) 33 { 34 switch (i) { 35 case 0: 36 switch (j) { 37 case 0: return 0; 38 case 1: 39 switch (k) { 40 case 0: return 0; 41 case 1: return 0; 42 case 2: return 1; 43 } 44 case 2: 45 switch (k) { 46 case 0: return 0; 47 case 1: return -1; 48 case 2: return 0; 49 } 50 } 51 case 1: 52 switch (j) { 53 case 0: 54 switch (k) { 55 case 0: return 0; 56 case 1: return 0; 57 case 2: return -1; 58 } 59 case 1: return 0; 60 case 2: 61 switch (k) { 62 case 0: return 1; 63 case 1: return 0; 64 case 2: return 0; 65 } 66 } 67 case 2: 68 switch (j) { 69 case 0: 70 switch (k) { 71 case 0: return 0; 72 case 1: return 1; 73 case 2: return 0; 74 } 75 case 1: 76 switch (k) { 77 case 0: return -1; 78 case 1: return 0; 79 case 2: return 0; 80 } 81 case 2: return 0; 82 } 83 } 84 return 0; 85 } 86 87 #undef __FUNCT__ 88 #define __FUNCT__ "DMPlexProjectRigidBody" 89 static PetscErrorCode DMPlexProjectRigidBody(PetscInt dim, PetscReal t, const PetscReal X[], PetscInt Nf, PetscScalar *mode, void *ctx) 90 { 91 PetscInt *ctxInt = (PetscInt *) ctx; 92 PetscInt dim2 = ctxInt[0]; 93 PetscInt d = ctxInt[1]; 94 PetscInt i, j, k = dim > 2 ? d - dim : d; 95 96 PetscFunctionBegin; 97 if (dim != dim2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input dimension %d does not match context dimension %d", dim, dim2); 98 for (i = 0; i < dim; i++) mode[i] = 0.; 99 if (d < dim) { 100 mode[d] = 1.; 101 } else { 102 for (i = 0; i < dim; i++) { 103 for (j = 0; j < dim; j++) { 104 mode[j] += epsilon(i, j, k)*X[i]; 105 } 106 } 107 } 108 PetscFunctionReturn(0); 109 } 110 111 #undef __FUNCT__ 112 #define __FUNCT__ "DMPlexCreateRigidBody" 113 /*@C 114 DMPlexCreateRigidBody - for the default global section, create rigid body modes from coordinates 115 116 Collective on DM 117 118 Input Arguments: 119 . dm - the DM 120 121 Output Argument: 122 . sp - the null space 123 124 Note: This is necessary to take account of Dirichlet conditions on the displacements 125 126 Level: advanced 127 128 .seealso: MatNullSpaceCreate() 129 @*/ 130 PetscErrorCode DMPlexCreateRigidBody(DM dm, MatNullSpace *sp) 131 { 132 MPI_Comm comm; 133 Vec mode[6]; 134 PetscSection section, globalSection; 135 PetscInt dim, dimEmbed, n, m, d, i, j; 136 PetscErrorCode ierr; 137 138 PetscFunctionBegin; 139 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 140 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 141 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 142 if (dim == 1) { 143 ierr = MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);CHKERRQ(ierr); 144 PetscFunctionReturn(0); 145 } 146 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 147 ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr); 148 ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr); 149 m = (dim*(dim+1))/2; 150 ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr); 151 ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr); 152 ierr = VecSetUp(mode[0]);CHKERRQ(ierr); 153 for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);} 154 for (d = 0; d < m; d++) { 155 PetscInt ctx[2]; 156 void *voidctx = (void *) (&ctx[0]); 157 PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody; 158 159 ctx[0] = dimEmbed; 160 ctx[1] = d; 161 ierr = DMProjectFunction(dm, 0.0, &func, &voidctx, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 162 } 163 for (i = 0; i < dim; ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);} 164 /* Orthonormalize system */ 165 for (i = dim; i < m; ++i) { 166 PetscScalar dots[6]; 167 168 ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr); 169 for (j = 0; j < i; ++j) dots[j] *= -1.0; 170 ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr); 171 ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr); 172 } 173 ierr = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);CHKERRQ(ierr); 174 for (i = 0; i< m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);} 175 PetscFunctionReturn(0); 176 } 177 178 #undef __FUNCT__ 179 #define __FUNCT__ "DMPlexSetMaxProjectionHeight" 180 /*@ 181 DMPlexSetMaxProjectionHeight - In DMPlexProjectXXXLocal() functions, the projected values of a basis function's dofs 182 are computed by associating the basis function with one of the mesh points in its transitively-closed support, and 183 evaluating the dual space basis of that point. A basis function is associated with the point in its 184 transitively-closed support whose mesh height is highest (w.r.t. DAG height), but not greater than the maximum 185 projection height, which is set with this function. By default, the maximum projection height is zero, which means 186 that only mesh cells are used to project basis functions. A height of one, for example, evaluates a cell-interior 187 basis functions using its cells dual space basis, but all other basis functions with the dual space basis of a face. 188 189 Input Parameters: 190 + dm - the DMPlex object 191 - height - the maximum projection height >= 0 192 193 Level: advanced 194 195 .seealso: DMPlexGetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal() 196 @*/ 197 PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height) 198 { 199 DM_Plex *plex = (DM_Plex *) dm->data; 200 201 PetscFunctionBegin; 202 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 203 plex->maxProjectionHeight = height; 204 PetscFunctionReturn(0); 205 } 206 207 #undef __FUNCT__ 208 #define __FUNCT__ "DMPlexGetMaxProjectionHeight" 209 /*@ 210 DMPlexGetMaxProjectionHeight - Get the maximum height (w.r.t. DAG) of mesh points used to evaluate dual bases in 211 DMPlexProjectXXXLocal() functions. 212 213 Input Parameters: 214 . dm - the DMPlex object 215 216 Output Parameters: 217 . height - the maximum projection height 218 219 Level: intermediate 220 221 .seealso: DMPlexSetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal() 222 @*/ 223 PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height) 224 { 225 DM_Plex *plex = (DM_Plex *) dm->data; 226 227 PetscFunctionBegin; 228 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 229 *height = plex->maxProjectionHeight; 230 PetscFunctionReturn(0); 231 } 232 233 #undef __FUNCT__ 234 #define __FUNCT__ "DMProjectFunctionLabelLocal_Plex" 235 PetscErrorCode DMProjectFunctionLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 236 { 237 PetscDualSpace *sp, *cellsp; 238 PetscInt *numComp; 239 PetscSection section; 240 PetscScalar *values; 241 PetscBool *fieldActive; 242 PetscInt numFields, dim, dimEmbed, spDim, totDim = 0, numValues, pStart, pEnd, cStart, cEnd, cEndInterior, f, d, v, i, comp, maxHeight, h; 243 PetscErrorCode ierr; 244 245 PetscFunctionBegin; 246 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 247 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 248 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 249 if (cEnd <= cStart) PetscFunctionReturn(0); 250 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 251 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 252 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 253 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 254 ierr = PetscMalloc2(numFields,&sp,numFields,&numComp);CHKERRQ(ierr); 255 ierr = DMPlexGetMaxProjectionHeight(dm,&maxHeight);CHKERRQ(ierr); 256 if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"maximum projection height %d not in [0, %d)\n", maxHeight,dim);} 257 if (maxHeight > 0) {ierr = PetscMalloc1(numFields,&cellsp);CHKERRQ(ierr);} 258 else {cellsp = sp;} 259 for (h = 0; h <= maxHeight; h++) { 260 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 261 if (!h) {pStart = cStart; pEnd = cEnd;} 262 if (pEnd <= pStart) continue; 263 totDim = 0; 264 for (f = 0; f < numFields; ++f) { 265 PetscObject obj; 266 PetscClassId id; 267 268 ierr = DMGetField(dm, f, &obj);CHKERRQ(ierr); 269 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 270 if (id == PETSCFE_CLASSID) { 271 PetscFE fe = (PetscFE) obj; 272 273 ierr = PetscFEGetNumComponents(fe, &numComp[f]);CHKERRQ(ierr); 274 if (!h) { 275 ierr = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr); 276 sp[f] = cellsp[f]; 277 ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr); 278 } else { 279 ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);CHKERRQ(ierr); 280 if (!sp[f]) continue; 281 } 282 } else if (id == PETSCFV_CLASSID) { 283 PetscFV fv = (PetscFV) obj; 284 285 ierr = PetscFVGetNumComponents(fv, &numComp[f]);CHKERRQ(ierr); 286 ierr = PetscFVGetDualSpace(fv, &sp[f]);CHKERRQ(ierr); 287 ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr); 288 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 289 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 290 totDim += spDim*numComp[f]; 291 } 292 ierr = DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);CHKERRQ(ierr); 293 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim); 294 if (!totDim) continue; 295 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 296 ierr = DMGetWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);CHKERRQ(ierr); 297 for (f = 0; f < numFields; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE; 298 for (i = 0; i < numIds; ++i) { 299 IS pointIS; 300 const PetscInt *points; 301 PetscInt n, p; 302 303 ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr); 304 ierr = ISGetLocalSize(pointIS, &n);CHKERRQ(ierr); 305 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 306 for (p = 0; p < n; ++p) { 307 const PetscInt point = points[p]; 308 PetscFECellGeom geom; 309 310 if ((point < pStart) || (point >= pEnd)) continue; 311 ierr = DMPlexComputeCellGeometryFEM(dm, point, NULL, geom.v0, geom.J, NULL, &geom.detJ);CHKERRQ(ierr); 312 geom.dim = dim - h; 313 geom.dimEmbed = dimEmbed; 314 for (f = 0, v = 0; f < numFields; ++f) { 315 void * const ctx = ctxs ? ctxs[f] : NULL; 316 317 if (!sp[f]) continue; 318 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 319 for (d = 0; d < spDim; ++d) { 320 if (funcs[f]) { 321 ierr = PetscDualSpaceApply(sp[f], d, time, &geom, numComp[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr); 322 } else { 323 for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0; 324 } 325 v += numComp[f]; 326 } 327 } 328 ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, values, mode);CHKERRQ(ierr); 329 } 330 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 331 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 332 } 333 if (h) { 334 for (f = 0; f < numFields; ++f) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);} 335 } 336 } 337 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 338 ierr = DMRestoreWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);CHKERRQ(ierr); 339 for (f = 0; f < numFields; ++f) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);} 340 ierr = PetscFree2(sp, numComp);CHKERRQ(ierr); 341 if (maxHeight > 0) { 342 ierr = PetscFree(cellsp);CHKERRQ(ierr); 343 } 344 PetscFunctionReturn(0); 345 } 346 347 #undef __FUNCT__ 348 #define __FUNCT__ "DMProjectFunctionLocal_Plex" 349 PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 350 { 351 PetscDualSpace *sp, *cellsp; 352 PetscInt *numComp; 353 PetscSection section; 354 PetscScalar *values; 355 PetscInt Nf, dim, dimEmbed, spDim, totDim = 0, numValues, pStart, pEnd, p, cStart, cEnd, cEndInterior, f, d, v, comp, h, maxHeight; 356 PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE; 357 PetscErrorCode ierr; 358 359 PetscFunctionBegin; 360 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 361 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 362 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 363 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 364 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 365 ierr = PetscMalloc3(Nf, &isFE, Nf, &sp, Nf, &numComp);CHKERRQ(ierr); 366 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 367 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 368 ierr = DMPlexGetMaxProjectionHeight(dm,&maxHeight);CHKERRQ(ierr); 369 if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"maximum projection height %d not in [0, %d)\n", maxHeight,dim);} 370 if (maxHeight > 0) { 371 ierr = PetscMalloc1(Nf,&cellsp);CHKERRQ(ierr); 372 } 373 else { 374 cellsp = sp; 375 } 376 for (h = 0; h <= maxHeight; h++) { 377 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 378 if (!h) {pStart = cStart; pEnd = cEnd;} 379 totDim = 0; 380 for (f = 0; f < Nf; ++f) { 381 PetscObject obj; 382 PetscClassId id; 383 384 ierr = DMGetField(dm, f, &obj);CHKERRQ(ierr); 385 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 386 if (id == PETSCFE_CLASSID) { 387 PetscFE fe = (PetscFE) obj; 388 389 hasFE = PETSC_TRUE; 390 isFE[f] = PETSC_TRUE; 391 ierr = PetscFEGetNumComponents(fe, &numComp[f]);CHKERRQ(ierr); 392 if (!h) { 393 ierr = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr); 394 sp[f] = cellsp[f]; 395 ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr); 396 } else { 397 ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);CHKERRQ(ierr); 398 if (!sp[f]) { 399 continue; 400 } 401 } 402 } else if (id == PETSCFV_CLASSID) { 403 PetscFV fv = (PetscFV) obj; 404 405 hasFV = PETSC_TRUE; 406 isFE[f] = PETSC_FALSE; 407 if (!h) { 408 ierr = PetscFVGetNumComponents(fv, &numComp[f]);CHKERRQ(ierr); 409 ierr = PetscFVGetDualSpace(fv, &cellsp[f]);CHKERRQ(ierr); 410 ierr = PetscObjectReference((PetscObject) cellsp[f]);CHKERRQ(ierr); 411 sp[f] = cellsp[f]; 412 } else { 413 sp[f] = NULL; 414 continue; /* FV doesn't have fields that can be evaluated at faces, corners, etc. */ 415 } 416 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 417 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 418 totDim += spDim*numComp[f]; 419 } 420 if (pEnd <= pStart) continue; 421 ierr = DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);CHKERRQ(ierr); 422 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim); 423 if (!totDim) continue; 424 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 425 for (p = pStart; p < pEnd; ++p) { 426 PetscFECellGeom fegeom; 427 PetscFVCellGeom fvgeom; 428 429 if (hasFE) { 430 ierr = DMPlexComputeCellGeometryFEM(dm, p, NULL, fegeom.v0, fegeom.J, NULL, &fegeom.detJ);CHKERRQ(ierr); 431 fegeom.dim = dim - h; 432 fegeom.dimEmbed = dimEmbed; 433 } 434 if (hasFV) {ierr = DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);CHKERRQ(ierr);} 435 for (f = 0, v = 0; f < Nf; ++f) { 436 void * const ctx = ctxs ? ctxs[f] : NULL; 437 438 if (!sp[f]) continue; 439 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 440 for (d = 0; d < spDim; ++d) { 441 if (funcs[f]) { 442 if (isFE[f]) {ierr = PetscDualSpaceApply(sp[f], d, time, &fegeom, numComp[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr);} 443 else {ierr = PetscDualSpaceApplyFVM(sp[f], d, time, &fvgeom, numComp[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr);} 444 } else { 445 for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0; 446 } 447 v += numComp[f]; 448 } 449 } 450 ierr = DMPlexVecSetClosure(dm, section, localX, p, values, mode);CHKERRQ(ierr); 451 } 452 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 453 if (h || !maxHeight) { 454 for (f = 0; f < Nf; f++) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);} 455 } 456 } 457 ierr = PetscFree3(isFE, sp, numComp);CHKERRQ(ierr); 458 if (maxHeight > 0) { 459 for (f = 0; f < Nf; f++) {ierr = PetscDualSpaceDestroy(&cellsp[f]);CHKERRQ(ierr);} 460 ierr = PetscFree(cellsp);CHKERRQ(ierr); 461 } 462 PetscFunctionReturn(0); 463 } 464 465 #undef __FUNCT__ 466 #define __FUNCT__ "DMProjectFieldLocal_Plex" 467 PetscErrorCode DMProjectFieldLocal_Plex(DM dm, Vec localU, 468 void (**funcs)(PetscInt, PetscInt, PetscInt, 469 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 470 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 471 PetscReal, const PetscReal[], PetscScalar[]), 472 InsertMode mode, Vec localX) 473 { 474 DM dmAux; 475 PetscDS prob, probAux = NULL; 476 Vec A; 477 PetscSection section, sectionAux = NULL; 478 PetscDualSpace *sp; 479 PetscInt *Ncf; 480 PetscScalar *values, *u, *u_x, *a, *a_x; 481 PetscReal *x, *v0, *J, *invJ, detJ; 482 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL; 483 PetscInt Nf, NfAux = 0, dim, spDim, totDim, numValues, cStart, cEnd, cEndInterior, c, f, d, v, comp, maxHeight; 484 PetscErrorCode ierr; 485 486 PetscFunctionBegin; 487 ierr = DMPlexGetMaxProjectionHeight(dm,&maxHeight);CHKERRQ(ierr); 488 if (maxHeight > 0) {SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Field projection for height > 0 not supported yet");} 489 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 490 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 491 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 492 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 493 ierr = PetscMalloc2(Nf, &sp, Nf, &Ncf);CHKERRQ(ierr); 494 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 495 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 496 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 497 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 498 ierr = PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);CHKERRQ(ierr); 499 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 500 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 501 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 502 if (dmAux) { 503 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 504 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 505 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 506 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 507 ierr = PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);CHKERRQ(ierr); 508 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);CHKERRQ(ierr); 509 } 510 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localU, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 511 ierr = DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr); 512 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim); 513 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 514 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 515 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 516 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 517 for (c = cStart; c < cEnd; ++c) { 518 PetscScalar *coefficients = NULL, *coefficientsAux = NULL; 519 520 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 521 ierr = DMPlexVecGetClosure(dm, section, localU, c, NULL, &coefficients);CHKERRQ(ierr); 522 if (dmAux) {ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);CHKERRQ(ierr);} 523 for (f = 0, v = 0; f < Nf; ++f) { 524 PetscObject obj; 525 PetscClassId id; 526 527 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 528 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 529 if (id == PETSCFE_CLASSID) { 530 PetscFE fe = (PetscFE) obj; 531 532 ierr = PetscFEGetDualSpace(fe, &sp[f]);CHKERRQ(ierr); 533 ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr); 534 ierr = PetscFEGetNumComponents(fe, &Ncf[f]);CHKERRQ(ierr); 535 } else if (id == PETSCFV_CLASSID) { 536 PetscFV fv = (PetscFV) obj; 537 538 ierr = PetscFVGetNumComponents(fv, &Ncf[f]);CHKERRQ(ierr); 539 ierr = PetscFVGetDualSpace(fv, &sp[f]);CHKERRQ(ierr); 540 ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr); 541 } 542 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 543 for (d = 0; d < spDim; ++d) { 544 PetscQuadrature quad; 545 const PetscReal *points, *weights; 546 PetscInt numPoints, q; 547 548 if (funcs[f]) { 549 ierr = PetscDualSpaceGetFunctional(sp[f], d, &quad);CHKERRQ(ierr); 550 ierr = PetscQuadratureGetData(quad, NULL, &numPoints, &points, &weights);CHKERRQ(ierr); 551 for (q = 0; q < numPoints; ++q) { 552 CoordinatesRefToReal(dim, dim, v0, J, &points[q*dim], x); 553 ierr = EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, coefficients, NULL, u, u_x, NULL);CHKERRQ(ierr); 554 ierr = EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, coefficientsAux, NULL, a, a_x, NULL);CHKERRQ(ierr); 555 (*funcs[f])(dim, Nf, NfAux, uOff, uOff_x, u, NULL, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, x, &values[v]); 556 } 557 } else { 558 for (comp = 0; comp < Ncf[f]; ++comp) values[v+comp] = 0.0; 559 } 560 v += Ncf[f]; 561 } 562 } 563 ierr = DMPlexVecRestoreClosure(dm, section, localU, c, NULL, &coefficients);CHKERRQ(ierr); 564 if (dmAux) {ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);CHKERRQ(ierr);} 565 ierr = DMPlexVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr); 566 } 567 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 568 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 569 for (f = 0; f < Nf; f++) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);} 570 ierr = PetscFree2(sp, Ncf);CHKERRQ(ierr); 571 PetscFunctionReturn(0); 572 } 573 574 #undef __FUNCT__ 575 #define __FUNCT__ "DMPlexInsertBoundaryValues_FEM_Internal" 576 static PetscErrorCode DMPlexInsertBoundaryValues_FEM_Internal(DM dm, PetscReal time, PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void *ctx, Vec locX) 577 { 578 PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx); 579 void **ctxs; 580 PetscInt numFields; 581 PetscErrorCode ierr; 582 583 PetscFunctionBegin; 584 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 585 ierr = PetscCalloc2(numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 586 funcs[field] = func; 587 ctxs[field] = ctx; 588 ierr = DMProjectFunctionLabelLocal(dm, time, label, numids, ids, funcs, ctxs, INSERT_BC_VALUES, locX);CHKERRQ(ierr); 589 ierr = PetscFree2(funcs,ctxs);CHKERRQ(ierr); 590 PetscFunctionReturn(0); 591 } 592 593 #undef __FUNCT__ 594 #define __FUNCT__ "DMPlexInsertBoundaryValues_FVM_Internal" 595 /* This ignores numcomps/comps */ 596 static PetscErrorCode DMPlexInsertBoundaryValues_FVM_Internal(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad, 597 PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX) 598 { 599 PetscDS prob; 600 PetscSF sf; 601 DM dmFace, dmCell, dmGrad; 602 const PetscScalar *facegeom, *cellgeom, *grad; 603 const PetscInt *leaves; 604 PetscScalar *x, *fx; 605 PetscInt dim, nleaves, loc, fStart, fEnd, pdim, i; 606 PetscErrorCode ierr; 607 608 PetscFunctionBegin; 609 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 610 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 611 nleaves = PetscMax(0, nleaves); 612 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 613 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 614 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 615 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 616 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 617 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 618 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 619 if (Grad) { 620 PetscFV fv; 621 622 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);CHKERRQ(ierr); 623 ierr = VecGetDM(Grad, &dmGrad);CHKERRQ(ierr); 624 ierr = VecGetArrayRead(Grad, &grad);CHKERRQ(ierr); 625 ierr = PetscFVGetNumComponents(fv, &pdim);CHKERRQ(ierr); 626 ierr = DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 627 } 628 ierr = VecGetArray(locX, &x);CHKERRQ(ierr); 629 for (i = 0; i < numids; ++i) { 630 IS faceIS; 631 const PetscInt *faces; 632 PetscInt numFaces, f; 633 634 ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr); 635 if (!faceIS) continue; /* No points with that id on this process */ 636 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 637 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 638 for (f = 0; f < numFaces; ++f) { 639 const PetscInt face = faces[f], *cells; 640 const PetscFVFaceGeom *fg; 641 642 if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */ 643 ierr = PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);CHKERRQ(ierr); 644 if (loc >= 0) continue; 645 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 646 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 647 if (Grad) { 648 const PetscFVCellGeom *cg; 649 const PetscScalar *cx, *cgrad; 650 PetscScalar *xG; 651 PetscReal dx[3]; 652 PetscInt d; 653 654 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);CHKERRQ(ierr); 655 ierr = DMPlexPointLocalRead(dm, cells[0], x, &cx);CHKERRQ(ierr); 656 ierr = DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);CHKERRQ(ierr); 657 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 658 DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx); 659 for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx); 660 ierr = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);CHKERRQ(ierr); 661 } else { 662 const PetscScalar *xI; 663 PetscScalar *xG; 664 665 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr); 666 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 667 ierr = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);CHKERRQ(ierr); 668 } 669 } 670 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 671 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 672 } 673 ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr); 674 if (Grad) { 675 ierr = DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 676 ierr = VecRestoreArrayRead(Grad, &grad);CHKERRQ(ierr); 677 } 678 ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 679 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 680 PetscFunctionReturn(0); 681 } 682 683 #undef __FUNCT__ 684 #define __FUNCT__ "DMPlexInsertBoundaryValues" 685 PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM) 686 { 687 PetscInt numBd, b; 688 PetscErrorCode ierr; 689 690 PetscFunctionBegin; 691 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 692 PetscValidHeaderSpecific(locX, VEC_CLASSID, 2); 693 if (faceGeomFVM) {PetscValidHeaderSpecific(faceGeomFVM, VEC_CLASSID, 4);} 694 if (cellGeomFVM) {PetscValidHeaderSpecific(cellGeomFVM, VEC_CLASSID, 5);} 695 if (gradFVM) {PetscValidHeaderSpecific(gradFVM, VEC_CLASSID, 6);} 696 ierr = DMGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 697 for (b = 0; b < numBd; ++b) { 698 PetscBool isEssential; 699 const char *labelname; 700 DMLabel label; 701 PetscInt field; 702 PetscObject obj; 703 PetscClassId id; 704 void (*func)(); 705 PetscInt numids; 706 const PetscInt *ids; 707 void *ctx; 708 709 ierr = DMGetBoundary(dm, b, &isEssential, NULL, &labelname, &field, NULL, NULL, &func, &numids, &ids, &ctx);CHKERRQ(ierr); 710 if (insertEssential != isEssential) continue; 711 ierr = DMGetLabel(dm, labelname, &label);CHKERRQ(ierr); 712 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 713 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 714 if (id == PETSCFE_CLASSID) { 715 if (!isEssential) continue; /* for FEM, there is no insertion to be done for non-essential boundary conditions */ 716 ierr = DMPlexInsertBoundaryValues_FEM_Internal(dm, time, field, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);CHKERRQ(ierr); 717 } else if (id == PETSCFV_CLASSID) { 718 if (!faceGeomFVM) continue; 719 ierr = DMPlexInsertBoundaryValues_FVM_Internal(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, 720 field, label, numids, ids, (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);CHKERRQ(ierr); 721 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 722 } 723 PetscFunctionReturn(0); 724 } 725 726 #undef __FUNCT__ 727 #define __FUNCT__ "DMComputeL2Diff_Plex" 728 PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 729 { 730 const PetscInt debug = 0; 731 PetscSection section; 732 PetscQuadrature quad; 733 Vec localX; 734 PetscScalar *funcVal, *interpolant; 735 PetscReal *coords, *v0, *J, *invJ, detJ; 736 PetscReal localDiff = 0.0; 737 const PetscReal *quadPoints, *quadWeights; 738 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 739 PetscErrorCode ierr; 740 741 PetscFunctionBegin; 742 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 743 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 744 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 745 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 746 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 747 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 748 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 749 for (field = 0; field < numFields; ++field) { 750 PetscObject obj; 751 PetscClassId id; 752 PetscInt Nc; 753 754 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 755 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 756 if (id == PETSCFE_CLASSID) { 757 PetscFE fe = (PetscFE) obj; 758 759 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 760 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 761 } else if (id == PETSCFV_CLASSID) { 762 PetscFV fv = (PetscFV) obj; 763 764 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 765 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 766 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 767 numComponents += Nc; 768 } 769 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 770 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 771 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 772 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 773 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 774 for (c = cStart; c < cEnd; ++c) { 775 PetscScalar *x = NULL; 776 PetscReal elemDiff = 0.0; 777 778 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 779 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 780 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 781 782 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 783 PetscObject obj; 784 PetscClassId id; 785 void * const ctx = ctxs ? ctxs[field] : NULL; 786 PetscInt Nb, Nc, q, fc; 787 788 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 789 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 790 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 791 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 792 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 793 if (debug) { 794 char title[1024]; 795 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 796 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 797 } 798 for (q = 0; q < numQuadPoints; ++q) { 799 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 800 ierr = (*funcs[field])(dim, time, coords, Nc, funcVal, ctx);CHKERRQ(ierr); 801 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 802 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 803 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 804 for (fc = 0; fc < Nc; ++fc) { 805 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ);CHKERRQ(ierr);} 806 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 807 } 808 } 809 fieldOffset += Nb*Nc; 810 } 811 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 812 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 813 localDiff += elemDiff; 814 } 815 ierr = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 816 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 817 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 818 *diff = PetscSqrtReal(*diff); 819 PetscFunctionReturn(0); 820 } 821 822 #undef __FUNCT__ 823 #define __FUNCT__ "DMComputeL2GradientDiff_Plex" 824 PetscErrorCode DMComputeL2GradientDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff) 825 { 826 const PetscInt debug = 0; 827 PetscSection section; 828 PetscQuadrature quad; 829 Vec localX; 830 PetscScalar *funcVal, *interpolantVec; 831 PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ; 832 PetscReal localDiff = 0.0; 833 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset, comp; 834 PetscErrorCode ierr; 835 836 PetscFunctionBegin; 837 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 838 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 839 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 840 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 841 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 842 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 843 for (field = 0; field < numFields; ++field) { 844 PetscFE fe; 845 PetscInt Nc; 846 847 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 848 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 849 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 850 numComponents += Nc; 851 } 852 /* ierr = DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 853 ierr = PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr); 854 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 855 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 856 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 857 for (c = cStart; c < cEnd; ++c) { 858 PetscScalar *x = NULL; 859 PetscReal elemDiff = 0.0; 860 861 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 862 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 863 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 864 865 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 866 PetscFE fe; 867 void * const ctx = ctxs ? ctxs[field] : NULL; 868 const PetscReal *quadPoints, *quadWeights; 869 PetscReal *basisDer; 870 PetscInt numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g; 871 872 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 873 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 874 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 875 ierr = PetscFEGetNumComponents(fe, &Ncomp);CHKERRQ(ierr); 876 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 877 if (debug) { 878 char title[1024]; 879 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 880 ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr); 881 } 882 for (q = 0; q < numQuadPoints; ++q) { 883 for (d = 0; d < dim; d++) { 884 coords[d] = v0[d]; 885 for (e = 0; e < dim; e++) { 886 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 887 } 888 } 889 ierr = (*funcs[field])(dim, time, coords, n, numFields, funcVal, ctx);CHKERRQ(ierr); 890 for (fc = 0; fc < Ncomp; ++fc) { 891 PetscScalar interpolant = 0.0; 892 893 for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0; 894 for (f = 0; f < Nb; ++f) { 895 const PetscInt fidx = f*Ncomp+fc; 896 897 for (d = 0; d < dim; ++d) { 898 realSpaceDer[d] = 0.0; 899 for (g = 0; g < dim; ++g) { 900 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g]; 901 } 902 interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d]; 903 } 904 } 905 for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d]; 906 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);CHKERRQ(ierr);} 907 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 908 } 909 } 910 comp += Ncomp; 911 fieldOffset += Nb*Ncomp; 912 } 913 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 914 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 915 localDiff += elemDiff; 916 } 917 ierr = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr); 918 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 919 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 920 *diff = PetscSqrtReal(*diff); 921 PetscFunctionReturn(0); 922 } 923 924 #undef __FUNCT__ 925 #define __FUNCT__ "DMComputeL2FieldDiff_Plex" 926 PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[]) 927 { 928 const PetscInt debug = 0; 929 PetscSection section; 930 PetscQuadrature quad; 931 Vec localX; 932 PetscScalar *funcVal, *interpolant; 933 PetscReal *coords, *v0, *J, *invJ, detJ; 934 PetscReal *localDiff; 935 const PetscReal *quadPoints, *quadWeights; 936 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 937 PetscErrorCode ierr; 938 939 PetscFunctionBegin; 940 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 941 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 942 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 943 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 944 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 945 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 946 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 947 for (field = 0; field < numFields; ++field) { 948 PetscObject obj; 949 PetscClassId id; 950 PetscInt Nc; 951 952 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 953 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 954 if (id == PETSCFE_CLASSID) { 955 PetscFE fe = (PetscFE) obj; 956 957 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 958 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 959 } else if (id == PETSCFV_CLASSID) { 960 PetscFV fv = (PetscFV) obj; 961 962 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 963 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 964 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 965 numComponents += Nc; 966 } 967 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 968 ierr = PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 969 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 970 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 971 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 972 for (c = cStart; c < cEnd; ++c) { 973 PetscScalar *x = NULL; 974 975 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 976 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 977 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 978 979 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 980 PetscObject obj; 981 PetscClassId id; 982 void * const ctx = ctxs ? ctxs[field] : NULL; 983 PetscInt Nb, Nc, q, fc; 984 985 PetscReal elemDiff = 0.0; 986 987 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 988 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 989 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 990 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 991 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 992 if (debug) { 993 char title[1024]; 994 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 995 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 996 } 997 for (q = 0; q < numQuadPoints; ++q) { 998 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 999 ierr = (*funcs[field])(dim, time, coords, numFields, funcVal, ctx);CHKERRQ(ierr); 1000 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1001 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1002 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1003 for (fc = 0; fc < Nc; ++fc) { 1004 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ);CHKERRQ(ierr);} 1005 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 1006 } 1007 } 1008 fieldOffset += Nb*Nc; 1009 localDiff[field] += elemDiff; 1010 } 1011 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1012 } 1013 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1014 ierr = MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1015 for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]); 1016 ierr = PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 1017 PetscFunctionReturn(0); 1018 } 1019 1020 #undef __FUNCT__ 1021 #define __FUNCT__ "DMPlexComputeL2DiffVec" 1022 /*@C 1023 DMPlexComputeL2DiffVec - This function computes the cellwise L_2 difference between a function u and an FEM interpolant solution u_h, and stores it in a Vec. 1024 1025 Input Parameters: 1026 + dm - The DM 1027 . time - The time 1028 . funcs - The functions to evaluate for each field component 1029 . ctxs - Optional array of contexts to pass to each function, or NULL. 1030 - X - The coefficient vector u_h 1031 1032 Output Parameter: 1033 . D - A Vec which holds the difference ||u - u_h||_2 for each cell 1034 1035 Level: developer 1036 1037 .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff() 1038 @*/ 1039 PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D) 1040 { 1041 PetscSection section; 1042 PetscQuadrature quad; 1043 Vec localX; 1044 PetscScalar *funcVal, *interpolant; 1045 PetscReal *coords, *v0, *J, *invJ, detJ; 1046 const PetscReal *quadPoints, *quadWeights; 1047 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 1048 PetscErrorCode ierr; 1049 1050 PetscFunctionBegin; 1051 ierr = VecSet(D, 0.0);CHKERRQ(ierr); 1052 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1053 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1054 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1055 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1056 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 1057 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1058 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1059 for (field = 0; field < numFields; ++field) { 1060 PetscObject obj; 1061 PetscClassId id; 1062 PetscInt Nc; 1063 1064 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1065 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1066 if (id == PETSCFE_CLASSID) { 1067 PetscFE fe = (PetscFE) obj; 1068 1069 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1070 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1071 } else if (id == PETSCFV_CLASSID) { 1072 PetscFV fv = (PetscFV) obj; 1073 1074 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1075 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1076 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1077 numComponents += Nc; 1078 } 1079 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 1080 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1081 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1082 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1083 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1084 for (c = cStart; c < cEnd; ++c) { 1085 PetscScalar *x = NULL; 1086 PetscScalar elemDiff = 0.0; 1087 1088 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1089 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 1090 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1091 1092 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1093 PetscObject obj; 1094 PetscClassId id; 1095 void * const ctx = ctxs ? ctxs[field] : NULL; 1096 PetscInt Nb, Nc, q, fc; 1097 1098 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1099 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1100 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1101 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1102 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1103 for (q = 0; q < numQuadPoints; ++q) { 1104 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 1105 ierr = (*funcs[field])(dim, time, coords, Nc, funcVal, ctx);CHKERRQ(ierr); 1106 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1107 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1108 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1109 for (fc = 0; fc < Nc; ++fc) { 1110 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 1111 } 1112 } 1113 fieldOffset += Nb*Nc; 1114 } 1115 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1116 ierr = VecSetValuesSection(D, section, c, &elemDiff, ADD_VALUES);CHKERRQ(ierr); 1117 } 1118 ierr = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 1119 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1120 ierr = VecSqrtAbs(D);CHKERRQ(ierr); 1121 PetscFunctionReturn(0); 1122 } 1123 1124 #undef __FUNCT__ 1125 #define __FUNCT__ "DMPlexComputeIntegralFEM" 1126 /*@ 1127 DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user 1128 1129 Input Parameters: 1130 + dm - The mesh 1131 . X - Local input vector 1132 - user - The user context 1133 1134 Output Parameter: 1135 . integral - Local integral for each field 1136 1137 Level: developer 1138 1139 .seealso: DMPlexComputeResidualFEM() 1140 @*/ 1141 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user) 1142 { 1143 DM_Plex *mesh = (DM_Plex *) dm->data; 1144 DM dmAux; 1145 Vec localX, A; 1146 PetscDS prob, probAux = NULL; 1147 PetscSection section, sectionAux; 1148 PetscFECellGeom *cgeom; 1149 PetscScalar *u, *a = NULL; 1150 PetscReal *lintegral, *vol; 1151 PetscInt dim, Nf, f, numCells, cStart, cEnd, cEndInterior, c; 1152 PetscInt totDim, totDimAux; 1153 PetscErrorCode ierr; 1154 1155 PetscFunctionBegin; 1156 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1157 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1158 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 1159 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1160 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1161 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1162 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1163 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1164 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1165 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1166 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1167 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1168 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1169 numCells = cEnd - cStart; 1170 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1171 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1172 if (dmAux) { 1173 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1174 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1175 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1176 } 1177 ierr = PetscMalloc4(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeom,numCells,&vol);CHKERRQ(ierr); 1178 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1179 for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;} 1180 for (c = cStart; c < cEnd; ++c) { 1181 PetscScalar *x = NULL; 1182 PetscInt i; 1183 1184 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeom[c].v0, cgeom[c].J, cgeom[c].invJ, &cgeom[c].detJ);CHKERRQ(ierr); 1185 ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol[c], NULL, NULL);CHKERRQ(ierr); 1186 if (cgeom[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeom[c].detJ, c); 1187 ierr = DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1188 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 1189 ierr = DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1190 if (dmAux) { 1191 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1192 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 1193 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1194 } 1195 } 1196 for (f = 0; f < Nf; ++f) { 1197 PetscObject obj; 1198 PetscClassId id; 1199 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 1200 1201 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1202 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1203 if (id == PETSCFE_CLASSID) { 1204 PetscFE fe = (PetscFE) obj; 1205 PetscQuadrature q; 1206 PetscInt Nq, Nb; 1207 1208 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1209 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1210 ierr = PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 1211 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1212 blockSize = Nb*Nq; 1213 batchSize = numBlocks * blockSize; 1214 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1215 numChunks = numCells / (numBatches*batchSize); 1216 Ne = numChunks*numBatches*batchSize; 1217 Nr = numCells % (numBatches*batchSize); 1218 offset = numCells - Nr; 1219 ierr = PetscFEIntegrate(fe, prob, f, Ne, cgeom, u, probAux, a, lintegral);CHKERRQ(ierr); 1220 ierr = PetscFEIntegrate(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);CHKERRQ(ierr); 1221 } else if (id == PETSCFV_CLASSID) { 1222 /* PetscFV fv = (PetscFV) obj; */ 1223 PetscInt foff; 1224 PetscPointFunc obj_func; 1225 PetscScalar lint; 1226 1227 ierr = PetscDSGetObjective(prob, f, &obj_func);CHKERRQ(ierr); 1228 ierr = PetscDSGetFieldOffset(prob, f, &foff);CHKERRQ(ierr); 1229 if (obj_func) { 1230 for (c = 0; c < numCells; ++c) { 1231 /* TODO: Need full pointwise interpolation and get centroid for x argument */ 1232 obj_func(dim, Nf, 0, NULL, NULL, &u[totDim*c+foff], NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0.0, NULL, &lint); 1233 lintegral[f] = PetscRealPart(lint)*vol[c]; 1234 } 1235 } 1236 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 1237 } 1238 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1239 if (mesh->printFEM) { 1240 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");CHKERRQ(ierr); 1241 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);CHKERRQ(ierr);} 1242 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 1243 } 1244 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1245 ierr = MPIU_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr); 1246 ierr = PetscFree4(lintegral,u,cgeom,vol);CHKERRQ(ierr); 1247 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1248 PetscFunctionReturn(0); 1249 } 1250 1251 #undef __FUNCT__ 1252 #define __FUNCT__ "DMPlexComputeInterpolatorNested" 1253 /*@ 1254 DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1255 1256 Input Parameters: 1257 + dmf - The fine mesh 1258 . dmc - The coarse mesh 1259 - user - The user context 1260 1261 Output Parameter: 1262 . In - The interpolation matrix 1263 1264 Level: developer 1265 1266 .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1267 @*/ 1268 PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user) 1269 { 1270 DM_Plex *mesh = (DM_Plex *) dmc->data; 1271 const char *name = "Interpolator"; 1272 PetscDS prob; 1273 PetscFE *feRef; 1274 PetscFV *fvRef; 1275 PetscSection fsection, fglobalSection; 1276 PetscSection csection, cglobalSection; 1277 PetscScalar *elemMat; 1278 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c; 1279 PetscInt cTotDim, rTotDim = 0; 1280 PetscErrorCode ierr; 1281 1282 PetscFunctionBegin; 1283 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1284 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1285 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1286 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1287 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1288 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1289 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1290 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1291 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1292 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1293 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr); 1294 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1295 for (f = 0; f < Nf; ++f) { 1296 PetscObject obj; 1297 PetscClassId id; 1298 PetscInt rNb = 0, Nc = 0; 1299 1300 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1301 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1302 if (id == PETSCFE_CLASSID) { 1303 PetscFE fe = (PetscFE) obj; 1304 1305 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1306 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1307 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1308 } else if (id == PETSCFV_CLASSID) { 1309 PetscFV fv = (PetscFV) obj; 1310 PetscDualSpace Q; 1311 1312 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1313 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1314 ierr = PetscDualSpaceGetDimension(Q, &rNb);CHKERRQ(ierr); 1315 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1316 } 1317 rTotDim += rNb*Nc; 1318 } 1319 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1320 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 1321 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 1322 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1323 PetscDualSpace Qref; 1324 PetscQuadrature f; 1325 const PetscReal *qpoints, *qweights; 1326 PetscReal *points; 1327 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1328 1329 /* Compose points from all dual basis functionals */ 1330 if (feRef[fieldI]) { 1331 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1332 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 1333 } else { 1334 ierr = PetscFVGetDualSpace(fvRef[fieldI], &Qref);CHKERRQ(ierr); 1335 ierr = PetscFVGetNumComponents(fvRef[fieldI], &Nc);CHKERRQ(ierr); 1336 } 1337 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1338 for (i = 0; i < fpdim; ++i) { 1339 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1340 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1341 npoints += Np; 1342 } 1343 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1344 for (i = 0, k = 0; i < fpdim; ++i) { 1345 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1346 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1347 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1348 } 1349 1350 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1351 PetscObject obj; 1352 PetscClassId id; 1353 PetscReal *B; 1354 PetscInt NcJ = 0, cpdim = 0, j; 1355 1356 ierr = PetscDSGetDiscretization(prob, fieldJ, &obj);CHKERRQ(ierr); 1357 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1358 if (id == PETSCFE_CLASSID) { 1359 PetscFE fe = (PetscFE) obj; 1360 1361 /* Evaluate basis at points */ 1362 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 1363 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1364 /* For now, fields only interpolate themselves */ 1365 if (fieldI == fieldJ) { 1366 if (Nc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in fine space field %d does not match coarse field %d", Nc, NcJ); 1367 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1368 for (i = 0, k = 0; i < fpdim; ++i) { 1369 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1370 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1371 for (p = 0; p < Np; ++p, ++k) { 1372 for (j = 0; j < cpdim; ++j) { 1373 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p]; 1374 } 1375 } 1376 } 1377 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1378 } 1379 } else if (id == PETSCFV_CLASSID) { 1380 PetscFV fv = (PetscFV) obj; 1381 1382 /* Evaluate constant function at points */ 1383 ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr); 1384 cpdim = 1; 1385 /* For now, fields only interpolate themselves */ 1386 if (fieldI == fieldJ) { 1387 if (Nc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in fine space field %d does not match coarse field %d", Nc, NcJ); 1388 for (i = 0, k = 0; i < fpdim; ++i) { 1389 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1390 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1391 for (p = 0; p < Np; ++p, ++k) { 1392 for (j = 0; j < cpdim; ++j) { 1393 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p]; 1394 } 1395 } 1396 } 1397 } 1398 } 1399 offsetJ += cpdim*NcJ; 1400 } 1401 offsetI += fpdim*Nc; 1402 ierr = PetscFree(points);CHKERRQ(ierr); 1403 } 1404 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1405 /* Preallocate matrix */ 1406 { 1407 Mat preallocator; 1408 PetscScalar *vals; 1409 PetscInt *cellCIndices, *cellFIndices; 1410 PetscInt locRows, locCols, cell; 1411 1412 ierr = MatGetLocalSize(In, &locRows, &locCols);CHKERRQ(ierr); 1413 ierr = MatCreate(PetscObjectComm((PetscObject) In), &preallocator);CHKERRQ(ierr); 1414 ierr = MatSetType(preallocator, MATPREALLOCATOR);CHKERRQ(ierr); 1415 ierr = MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1416 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 1417 ierr = PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1418 for (cell = cStart; cell < cEnd; ++cell) { 1419 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1420 ierr = MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);CHKERRQ(ierr); 1421 } 1422 ierr = PetscFree3(vals,cellCIndices,cellFIndices);CHKERRQ(ierr); 1423 ierr = MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1424 ierr = MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1425 ierr = MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);CHKERRQ(ierr); 1426 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 1427 } 1428 /* Fill matrix */ 1429 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1430 for (c = cStart; c < cEnd; ++c) { 1431 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1432 } 1433 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1434 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1435 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1436 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1437 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1438 if (mesh->printFEM) { 1439 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1440 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1441 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1442 } 1443 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1444 PetscFunctionReturn(0); 1445 } 1446 1447 #undef __FUNCT__ 1448 #define __FUNCT__ "DMPlexComputeInterpolatorGeneral" 1449 /*@ 1450 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 1451 1452 Input Parameters: 1453 + dmf - The fine mesh 1454 . dmc - The coarse mesh 1455 - user - The user context 1456 1457 Output Parameter: 1458 . In - The interpolation matrix 1459 1460 Level: developer 1461 1462 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1463 @*/ 1464 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 1465 { 1466 DM_Plex *mesh = (DM_Plex *) dmf->data; 1467 const char *name = "Interpolator"; 1468 PetscDS prob; 1469 PetscSection fsection, csection, globalFSection, globalCSection; 1470 PetscHashJK ht; 1471 PetscLayout rLayout; 1472 PetscInt *dnz, *onz; 1473 PetscInt locRows, rStart, rEnd; 1474 PetscReal *x, *v0, *J, *invJ, detJ; 1475 PetscReal *v0c, *Jc, *invJc, detJc; 1476 PetscScalar *elemMat; 1477 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1478 PetscErrorCode ierr; 1479 1480 PetscFunctionBegin; 1481 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1482 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1483 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1484 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1485 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1486 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1487 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1488 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1489 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1490 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1491 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1492 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1493 ierr = PetscMalloc1(totDim*totDim, &elemMat);CHKERRQ(ierr); 1494 1495 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1496 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1497 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1498 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1499 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1500 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1501 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1502 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1503 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1504 for (field = 0; field < Nf; ++field) { 1505 PetscObject obj; 1506 PetscClassId id; 1507 PetscDualSpace Q = NULL; 1508 PetscQuadrature f; 1509 const PetscReal *qpoints; 1510 PetscInt Nc, Np, fpdim, i, d; 1511 1512 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1513 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1514 if (id == PETSCFE_CLASSID) { 1515 PetscFE fe = (PetscFE) obj; 1516 1517 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1518 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1519 } else if (id == PETSCFV_CLASSID) { 1520 PetscFV fv = (PetscFV) obj; 1521 1522 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1523 Nc = 1; 1524 } 1525 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1526 /* For each fine grid cell */ 1527 for (cell = cStart; cell < cEnd; ++cell) { 1528 PetscInt *findices, *cindices; 1529 PetscInt numFIndices, numCIndices; 1530 1531 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1532 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1533 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1534 for (i = 0; i < fpdim; ++i) { 1535 Vec pointVec; 1536 PetscScalar *pV; 1537 IS coarseCellIS; 1538 const PetscInt *coarseCells; 1539 PetscInt numCoarseCells, q, r, c; 1540 1541 /* Get points from the dual basis functional quadrature */ 1542 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1543 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1544 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1545 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1546 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1547 for (q = 0; q < Np; ++q) { 1548 /* Transform point to real space */ 1549 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1550 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1551 } 1552 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1553 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1554 ierr = DMLocatePoints(dmc, pointVec, &coarseCellIS);CHKERRQ(ierr); 1555 ierr = ISViewFromOptions(coarseCellIS, NULL, "-interp_is_view");CHKERRQ(ierr); 1556 /* Update preallocation info */ 1557 ierr = ISGetLocalSize(coarseCellIS, &numCoarseCells);CHKERRQ(ierr); 1558 ierr = ISGetIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1559 for (r = 0; r < Nc; ++r) { 1560 PetscHashJKKey key; 1561 PetscHashJKIter missing, iter; 1562 1563 key.j = findices[i*Nc+r]; 1564 if (key.j < 0) continue; 1565 /* Get indices for coarse elements */ 1566 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1567 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1568 for (c = 0; c < numCIndices; ++c) { 1569 key.k = cindices[c]; 1570 if (key.k < 0) continue; 1571 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1572 if (missing) { 1573 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1574 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1575 else ++onz[key.j-rStart]; 1576 } 1577 } 1578 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1579 } 1580 } 1581 ierr = ISRestoreIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1582 ierr = ISDestroy(&coarseCellIS);CHKERRQ(ierr); 1583 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1584 } 1585 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1586 } 1587 } 1588 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1589 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1590 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1591 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1592 for (field = 0; field < Nf; ++field) { 1593 PetscObject obj; 1594 PetscClassId id; 1595 PetscDualSpace Q = NULL; 1596 PetscQuadrature f; 1597 const PetscReal *qpoints, *qweights; 1598 PetscInt Nc, Np, fpdim, i, d; 1599 1600 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1601 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1602 if (id == PETSCFE_CLASSID) { 1603 PetscFE fe = (PetscFE) obj; 1604 1605 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1606 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1607 } else if (id == PETSCFV_CLASSID) { 1608 PetscFV fv = (PetscFV) obj; 1609 1610 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1611 Nc = 1; 1612 } 1613 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1614 /* For each fine grid cell */ 1615 for (cell = cStart; cell < cEnd; ++cell) { 1616 PetscInt *findices, *cindices; 1617 PetscInt numFIndices, numCIndices; 1618 1619 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1620 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1621 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1622 for (i = 0; i < fpdim; ++i) { 1623 Vec pointVec; 1624 PetscScalar *pV; 1625 IS coarseCellIS; 1626 const PetscInt *coarseCells; 1627 PetscInt numCoarseCells, cpdim, q, c, j; 1628 1629 /* Get points from the dual basis functional quadrature */ 1630 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1631 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, &qweights);CHKERRQ(ierr); 1632 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1633 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1634 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1635 for (q = 0; q < Np; ++q) { 1636 /* Transform point to real space */ 1637 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1638 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1639 } 1640 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1641 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1642 ierr = DMLocatePoints(dmc, pointVec, &coarseCellIS);CHKERRQ(ierr); 1643 /* Update preallocation info */ 1644 ierr = ISGetLocalSize(coarseCellIS, &numCoarseCells);CHKERRQ(ierr); 1645 ierr = ISGetIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1646 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1647 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1648 PetscReal pVReal[3]; 1649 1650 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1651 /* Transform points from real space to coarse reference space */ 1652 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell], NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 1653 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 1654 CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x); 1655 1656 if (id == PETSCFE_CLASSID) { 1657 PetscFE fe = (PetscFE) obj; 1658 PetscReal *B; 1659 1660 /* Evaluate coarse basis on contained point */ 1661 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1662 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 1663 /* Get elemMat entries by multiplying by weight */ 1664 for (j = 0; j < cpdim; ++j) { 1665 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = B[j*Nc + c]*qweights[ccell]; 1666 } 1667 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1668 } else { 1669 cpdim = 1; 1670 for (j = 0; j < cpdim; ++j) { 1671 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = 1.0*qweights[ccell]; 1672 } 1673 } 1674 /* Update interpolator */ 1675 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, Nc, numCIndices, elemMat);CHKERRQ(ierr);} 1676 ierr = MatSetValues(In, Nc, &findices[i*Nc], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1677 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1678 } 1679 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1680 ierr = ISRestoreIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1681 ierr = ISDestroy(&coarseCellIS);CHKERRQ(ierr); 1682 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1683 } 1684 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1685 } 1686 } 1687 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 1688 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 1689 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1690 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1691 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1692 PetscFunctionReturn(0); 1693 } 1694 1695 #undef __FUNCT__ 1696 #define __FUNCT__ "DMPlexComputeInjectorFEM" 1697 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 1698 { 1699 PetscDS prob; 1700 PetscFE *feRef; 1701 PetscFV *fvRef; 1702 Vec fv, cv; 1703 IS fis, cis; 1704 PetscSection fsection, fglobalSection, csection, cglobalSection; 1705 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 1706 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m; 1707 PetscErrorCode ierr; 1708 1709 PetscFunctionBegin; 1710 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1711 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1712 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1713 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1714 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1715 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1716 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1717 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1718 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1719 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1720 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1721 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1722 for (f = 0; f < Nf; ++f) { 1723 PetscObject obj; 1724 PetscClassId id; 1725 PetscInt fNb = 0, Nc = 0; 1726 1727 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1728 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1729 if (id == PETSCFE_CLASSID) { 1730 PetscFE fe = (PetscFE) obj; 1731 1732 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1733 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 1734 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1735 } else if (id == PETSCFV_CLASSID) { 1736 PetscFV fv = (PetscFV) obj; 1737 PetscDualSpace Q; 1738 1739 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1740 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1741 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 1742 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1743 } 1744 fTotDim += fNb*Nc; 1745 } 1746 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1747 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 1748 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 1749 PetscFE feC; 1750 PetscFV fvC; 1751 PetscDualSpace QF, QC; 1752 PetscInt NcF, NcC, fpdim, cpdim; 1753 1754 if (feRef[field]) { 1755 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 1756 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 1757 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 1758 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 1759 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1760 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 1761 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1762 } else { 1763 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 1764 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 1765 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 1766 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 1767 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1768 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 1769 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1770 } 1771 if (NcF != NcC) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in fine space field %d does not match coarse field %d", NcF, NcC); 1772 for (c = 0; c < cpdim; ++c) { 1773 PetscQuadrature cfunc; 1774 const PetscReal *cqpoints; 1775 PetscInt NpC; 1776 PetscBool found = PETSC_FALSE; 1777 1778 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 1779 ierr = PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 1780 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 1781 for (f = 0; f < fpdim; ++f) { 1782 PetscQuadrature ffunc; 1783 const PetscReal *fqpoints; 1784 PetscReal sum = 0.0; 1785 PetscInt NpF, comp; 1786 1787 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 1788 ierr = PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 1789 if (NpC != NpF) continue; 1790 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 1791 if (sum > 1.0e-9) continue; 1792 for (comp = 0; comp < NcC; ++comp) { 1793 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 1794 } 1795 found = PETSC_TRUE; 1796 break; 1797 } 1798 if (!found) { 1799 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 1800 if (fvRef[field]) { 1801 PetscInt comp; 1802 for (comp = 0; comp < NcC; ++comp) { 1803 cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp; 1804 } 1805 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 1806 } 1807 } 1808 offsetC += cpdim*NcC; 1809 offsetF += fpdim*NcF; 1810 } 1811 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 1812 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1813 1814 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 1815 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 1816 ierr = VecGetOwnershipRange(cv, &startC, &endC);CHKERRQ(ierr); 1817 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 1818 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 1819 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 1820 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 1821 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 1822 for (c = cStart; c < cEnd; ++c) { 1823 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 1824 for (d = 0; d < cTotDim; ++d) { 1825 if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue; 1826 if ((findices[cellCIndices[d]-startC] >= 0) && (findices[cellCIndices[d]-startC] != cellFIndices[cmap[d]])) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Coarse dof %d maps to both %d and %d", cindices[cellCIndices[d]-startC], findices[cellCIndices[d]-startC], cellFIndices[cmap[d]]); 1827 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 1828 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 1829 } 1830 } 1831 ierr = PetscFree(cmap);CHKERRQ(ierr); 1832 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 1833 1834 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 1835 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 1836 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 1837 ierr = ISDestroy(&cis);CHKERRQ(ierr); 1838 ierr = ISDestroy(&fis);CHKERRQ(ierr); 1839 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 1840 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 1841 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1842 PetscFunctionReturn(0); 1843 } 1844