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" 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, 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, 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 (!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 ierr = DMPlexInsertBoundaryValues_FEM_Internal(dm, time, field, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);CHKERRQ(ierr); 716 } else if (id == PETSCFV_CLASSID) { 717 if (!faceGeomFVM) continue; 718 ierr = DMPlexInsertBoundaryValues_FVM_Internal(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, 719 field, label, numids, ids, (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);CHKERRQ(ierr); 720 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 721 } 722 PetscFunctionReturn(0); 723 } 724 725 #undef __FUNCT__ 726 #define __FUNCT__ "DMComputeL2Diff_Plex" 727 PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 728 { 729 const PetscInt debug = 0; 730 PetscSection section; 731 PetscQuadrature quad; 732 Vec localX; 733 PetscScalar *funcVal, *interpolant; 734 PetscReal *coords, *v0, *J, *invJ, detJ; 735 PetscReal localDiff = 0.0; 736 const PetscReal *quadPoints, *quadWeights; 737 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 738 PetscErrorCode ierr; 739 740 PetscFunctionBegin; 741 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 742 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 743 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 744 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 745 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 746 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 747 for (field = 0; field < numFields; ++field) { 748 PetscObject obj; 749 PetscClassId id; 750 PetscInt Nc; 751 752 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 753 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 754 if (id == PETSCFE_CLASSID) { 755 PetscFE fe = (PetscFE) obj; 756 757 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 758 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 759 } else if (id == PETSCFV_CLASSID) { 760 PetscFV fv = (PetscFV) obj; 761 762 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 763 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 764 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 765 numComponents += Nc; 766 } 767 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 768 ierr = DMProjectFunctionLocal(dm, 0.0, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 769 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 770 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 771 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 772 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 773 for (c = cStart; c < cEnd; ++c) { 774 PetscScalar *x = NULL; 775 PetscReal elemDiff = 0.0; 776 777 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 778 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 779 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 780 781 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 782 PetscObject obj; 783 PetscClassId id; 784 void * const ctx = ctxs ? ctxs[field] : NULL; 785 PetscInt Nb, Nc, q, fc; 786 787 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 788 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 789 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 790 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 791 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 792 if (debug) { 793 char title[1024]; 794 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 795 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 796 } 797 for (q = 0; q < numQuadPoints; ++q) { 798 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 799 ierr = (*funcs[field])(dim, time, coords, Nc, funcVal, ctx);CHKERRQ(ierr); 800 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 801 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 802 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 803 for (fc = 0; fc < Nc; ++fc) { 804 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);} 805 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 806 } 807 } 808 fieldOffset += Nb*Nc; 809 } 810 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 811 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 812 localDiff += elemDiff; 813 } 814 ierr = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 815 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 816 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 817 *diff = PetscSqrtReal(*diff); 818 PetscFunctionReturn(0); 819 } 820 821 #undef __FUNCT__ 822 #define __FUNCT__ "DMComputeL2GradientDiff_Plex" 823 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) 824 { 825 const PetscInt debug = 0; 826 PetscSection section; 827 PetscQuadrature quad; 828 Vec localX; 829 PetscScalar *funcVal, *interpolantVec; 830 PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ; 831 PetscReal localDiff = 0.0; 832 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset, comp; 833 PetscErrorCode ierr; 834 835 PetscFunctionBegin; 836 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 837 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 838 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 839 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 840 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 841 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 842 for (field = 0; field < numFields; ++field) { 843 PetscFE fe; 844 PetscInt Nc; 845 846 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 847 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 848 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 849 numComponents += Nc; 850 } 851 /* ierr = DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 852 ierr = PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr); 853 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 854 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 855 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 856 for (c = cStart; c < cEnd; ++c) { 857 PetscScalar *x = NULL; 858 PetscReal elemDiff = 0.0; 859 860 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 861 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 862 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 863 864 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 865 PetscFE fe; 866 void * const ctx = ctxs ? ctxs[field] : NULL; 867 const PetscReal *quadPoints, *quadWeights; 868 PetscReal *basisDer; 869 PetscInt numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g; 870 871 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 872 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 873 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 874 ierr = PetscFEGetNumComponents(fe, &Ncomp);CHKERRQ(ierr); 875 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 876 if (debug) { 877 char title[1024]; 878 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 879 ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr); 880 } 881 for (q = 0; q < numQuadPoints; ++q) { 882 for (d = 0; d < dim; d++) { 883 coords[d] = v0[d]; 884 for (e = 0; e < dim; e++) { 885 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 886 } 887 } 888 ierr = (*funcs[field])(dim, time, coords, n, numFields, funcVal, ctx);CHKERRQ(ierr); 889 for (fc = 0; fc < Ncomp; ++fc) { 890 PetscScalar interpolant = 0.0; 891 892 for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0; 893 for (f = 0; f < Nb; ++f) { 894 const PetscInt fidx = f*Ncomp+fc; 895 896 for (d = 0; d < dim; ++d) { 897 realSpaceDer[d] = 0.0; 898 for (g = 0; g < dim; ++g) { 899 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g]; 900 } 901 interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d]; 902 } 903 } 904 for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d]; 905 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);} 906 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 907 } 908 } 909 comp += Ncomp; 910 fieldOffset += Nb*Ncomp; 911 } 912 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 913 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 914 localDiff += elemDiff; 915 } 916 ierr = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr); 917 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 918 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 919 *diff = PetscSqrtReal(*diff); 920 PetscFunctionReturn(0); 921 } 922 923 #undef __FUNCT__ 924 #define __FUNCT__ "DMComputeL2FieldDiff_Plex" 925 PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[]) 926 { 927 const PetscInt debug = 0; 928 PetscSection section; 929 PetscQuadrature quad; 930 Vec localX; 931 PetscScalar *funcVal, *interpolant; 932 PetscReal *coords, *v0, *J, *invJ, detJ; 933 PetscReal *localDiff; 934 const PetscReal *quadPoints, *quadWeights; 935 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 936 PetscErrorCode ierr; 937 938 PetscFunctionBegin; 939 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 940 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 941 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 942 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 943 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 944 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 945 for (field = 0; field < numFields; ++field) { 946 PetscObject obj; 947 PetscClassId id; 948 PetscInt Nc; 949 950 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 951 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 952 if (id == PETSCFE_CLASSID) { 953 PetscFE fe = (PetscFE) obj; 954 955 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 956 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 957 } else if (id == PETSCFV_CLASSID) { 958 PetscFV fv = (PetscFV) obj; 959 960 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 961 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 962 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 963 numComponents += Nc; 964 } 965 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 966 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 967 ierr = PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 968 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 969 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 970 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 971 for (c = cStart; c < cEnd; ++c) { 972 PetscScalar *x = NULL; 973 974 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 975 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 976 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 977 978 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 979 PetscObject obj; 980 PetscClassId id; 981 void * const ctx = ctxs ? ctxs[field] : NULL; 982 PetscInt Nb, Nc, q, fc; 983 984 PetscReal elemDiff = 0.0; 985 986 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 987 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 988 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 989 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 990 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 991 if (debug) { 992 char title[1024]; 993 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 994 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 995 } 996 for (q = 0; q < numQuadPoints; ++q) { 997 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 998 ierr = (*funcs[field])(dim, time, coords, numFields, funcVal, ctx);CHKERRQ(ierr); 999 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1000 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1001 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1002 for (fc = 0; fc < Nc; ++fc) { 1003 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);} 1004 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 1005 } 1006 } 1007 fieldOffset += Nb*Nc; 1008 localDiff[field] += elemDiff; 1009 } 1010 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1011 } 1012 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1013 ierr = MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1014 for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]); 1015 ierr = PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 1016 PetscFunctionReturn(0); 1017 } 1018 1019 #undef __FUNCT__ 1020 #define __FUNCT__ "DMPlexComputeL2DiffVec" 1021 /*@C 1022 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. 1023 1024 Input Parameters: 1025 + dm - The DM 1026 . time - The time 1027 . funcs - The functions to evaluate for each field component 1028 . ctxs - Optional array of contexts to pass to each function, or NULL. 1029 - X - The coefficient vector u_h 1030 1031 Output Parameter: 1032 . D - A Vec which holds the difference ||u - u_h||_2 for each cell 1033 1034 Level: developer 1035 1036 .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff() 1037 @*/ 1038 PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D) 1039 { 1040 PetscSection section; 1041 PetscQuadrature quad; 1042 Vec localX; 1043 PetscScalar *funcVal, *interpolant; 1044 PetscReal *coords, *v0, *J, *invJ, detJ; 1045 const PetscReal *quadPoints, *quadWeights; 1046 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 1047 PetscErrorCode ierr; 1048 1049 PetscFunctionBegin; 1050 ierr = VecSet(D, 0.0);CHKERRQ(ierr); 1051 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1052 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1053 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1054 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1055 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1056 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1057 for (field = 0; field < numFields; ++field) { 1058 PetscObject obj; 1059 PetscClassId id; 1060 PetscInt Nc; 1061 1062 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1063 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1064 if (id == PETSCFE_CLASSID) { 1065 PetscFE fe = (PetscFE) obj; 1066 1067 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1068 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1069 } else if (id == PETSCFV_CLASSID) { 1070 PetscFV fv = (PetscFV) obj; 1071 1072 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1073 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1074 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1075 numComponents += Nc; 1076 } 1077 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 1078 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 1079 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1080 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1081 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1082 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1083 for (c = cStart; c < cEnd; ++c) { 1084 PetscScalar *x = NULL; 1085 PetscScalar elemDiff = 0.0; 1086 1087 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1088 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 1089 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1090 1091 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1092 PetscObject obj; 1093 PetscClassId id; 1094 void * const ctx = ctxs ? ctxs[field] : NULL; 1095 PetscInt Nb, Nc, q, fc; 1096 1097 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1098 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1099 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1100 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1101 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1102 for (q = 0; q < numQuadPoints; ++q) { 1103 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 1104 ierr = (*funcs[field])(dim, time, coords, Nc, funcVal, ctx);CHKERRQ(ierr); 1105 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1106 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1107 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1108 for (fc = 0; fc < Nc; ++fc) { 1109 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 1110 } 1111 } 1112 fieldOffset += Nb*Nc; 1113 } 1114 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1115 ierr = VecSetValuesSection(D, section, c, &elemDiff, ADD_VALUES);CHKERRQ(ierr); 1116 } 1117 ierr = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 1118 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1119 ierr = VecSqrtAbs(D);CHKERRQ(ierr); 1120 PetscFunctionReturn(0); 1121 } 1122 1123 #undef __FUNCT__ 1124 #define __FUNCT__ "DMPlexComputeIntegralFEM" 1125 /*@ 1126 DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user 1127 1128 Input Parameters: 1129 + dm - The mesh 1130 . X - Local input vector 1131 - user - The user context 1132 1133 Output Parameter: 1134 . integral - Local integral for each field 1135 1136 Level: developer 1137 1138 .seealso: DMPlexComputeResidualFEM() 1139 @*/ 1140 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user) 1141 { 1142 DM_Plex *mesh = (DM_Plex *) dm->data; 1143 DM dmAux; 1144 Vec localX, A; 1145 PetscDS prob, probAux = NULL; 1146 PetscSection section, sectionAux; 1147 PetscFECellGeom *cgeom; 1148 PetscScalar *u, *a = NULL; 1149 PetscReal *lintegral, *vol; 1150 PetscInt dim, Nf, f, numCells, cStart, cEnd, cEndInterior, c; 1151 PetscInt totDim, totDimAux; 1152 PetscErrorCode ierr; 1153 1154 PetscFunctionBegin; 1155 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1156 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1157 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1158 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1159 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1160 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1161 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1162 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1163 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1164 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1165 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1166 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1167 numCells = cEnd - cStart; 1168 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1169 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1170 if (dmAux) { 1171 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1172 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1173 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1174 } 1175 ierr = DMPlexInsertBoundaryValues(dm, localX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 1176 ierr = PetscMalloc4(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeom,numCells,&vol);CHKERRQ(ierr); 1177 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1178 for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;} 1179 for (c = cStart; c < cEnd; ++c) { 1180 PetscScalar *x = NULL; 1181 PetscInt i; 1182 1183 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeom[c].v0, cgeom[c].J, cgeom[c].invJ, &cgeom[c].detJ);CHKERRQ(ierr); 1184 ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol[c], NULL, NULL);CHKERRQ(ierr); 1185 if (cgeom[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeom[c].detJ, c); 1186 ierr = DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1187 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 1188 ierr = DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1189 if (dmAux) { 1190 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1191 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 1192 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1193 } 1194 } 1195 for (f = 0; f < Nf; ++f) { 1196 PetscObject obj; 1197 PetscClassId id; 1198 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 1199 1200 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1201 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1202 if (id == PETSCFE_CLASSID) { 1203 PetscFE fe = (PetscFE) obj; 1204 PetscQuadrature q; 1205 PetscInt Nq, Nb; 1206 1207 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1208 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1209 ierr = PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 1210 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1211 blockSize = Nb*Nq; 1212 batchSize = numBlocks * blockSize; 1213 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1214 numChunks = numCells / (numBatches*batchSize); 1215 Ne = numChunks*numBatches*batchSize; 1216 Nr = numCells % (numBatches*batchSize); 1217 offset = numCells - Nr; 1218 ierr = PetscFEIntegrate(fe, prob, f, Ne, cgeom, u, probAux, a, lintegral);CHKERRQ(ierr); 1219 ierr = PetscFEIntegrate(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);CHKERRQ(ierr); 1220 } else if (id == PETSCFV_CLASSID) { 1221 /* PetscFV fv = (PetscFV) obj; */ 1222 PetscInt foff; 1223 PetscPointFunc obj_func; 1224 PetscScalar lint; 1225 1226 ierr = PetscDSGetObjective(prob, f, &obj_func);CHKERRQ(ierr); 1227 ierr = PetscDSGetFieldOffset(prob, f, &foff);CHKERRQ(ierr); 1228 if (obj_func) { 1229 for (c = 0; c < numCells; ++c) { 1230 /* TODO: Need full pointwise interpolation and get centroid for x argument */ 1231 obj_func(dim, Nf, 0, NULL, NULL, &u[totDim*c+foff], NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0.0, NULL, &lint); 1232 lintegral[f] = PetscRealPart(lint)*vol[c]; 1233 } 1234 } 1235 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 1236 } 1237 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1238 if (mesh->printFEM) { 1239 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");CHKERRQ(ierr); 1240 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);CHKERRQ(ierr);} 1241 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 1242 } 1243 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1244 ierr = MPIU_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr); 1245 ierr = PetscFree4(lintegral,u,cgeom,vol);CHKERRQ(ierr); 1246 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1247 PetscFunctionReturn(0); 1248 } 1249 1250 #undef __FUNCT__ 1251 #define __FUNCT__ "DMPlexComputeInterpolatorNested" 1252 /*@ 1253 DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1254 1255 Input Parameters: 1256 + dmf - The fine mesh 1257 . dmc - The coarse mesh 1258 - user - The user context 1259 1260 Output Parameter: 1261 . In - The interpolation matrix 1262 1263 Level: developer 1264 1265 .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1266 @*/ 1267 PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user) 1268 { 1269 DM_Plex *mesh = (DM_Plex *) dmc->data; 1270 const char *name = "Interpolator"; 1271 PetscDS prob; 1272 PetscFE *feRef; 1273 PetscFV *fvRef; 1274 PetscSection fsection, fglobalSection; 1275 PetscSection csection, cglobalSection; 1276 PetscScalar *elemMat; 1277 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c; 1278 PetscInt cTotDim, rTotDim = 0; 1279 PetscErrorCode ierr; 1280 1281 PetscFunctionBegin; 1282 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1283 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1284 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1285 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1286 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1287 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1288 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1289 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1290 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1291 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1292 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr); 1293 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1294 for (f = 0; f < Nf; ++f) { 1295 PetscObject obj; 1296 PetscClassId id; 1297 PetscInt rNb = 0, Nc = 0; 1298 1299 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1300 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1301 if (id == PETSCFE_CLASSID) { 1302 PetscFE fe = (PetscFE) obj; 1303 1304 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1305 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1306 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1307 } else if (id == PETSCFV_CLASSID) { 1308 PetscFV fv = (PetscFV) obj; 1309 PetscDualSpace Q; 1310 1311 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1312 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1313 ierr = PetscDualSpaceGetDimension(Q, &rNb);CHKERRQ(ierr); 1314 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1315 } 1316 rTotDim += rNb*Nc; 1317 } 1318 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1319 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 1320 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 1321 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1322 PetscDualSpace Qref; 1323 PetscQuadrature f; 1324 const PetscReal *qpoints, *qweights; 1325 PetscReal *points; 1326 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1327 1328 /* Compose points from all dual basis functionals */ 1329 if (feRef[fieldI]) { 1330 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1331 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 1332 } else { 1333 ierr = PetscFVGetDualSpace(fvRef[fieldI], &Qref);CHKERRQ(ierr); 1334 ierr = PetscFVGetNumComponents(fvRef[fieldI], &Nc);CHKERRQ(ierr); 1335 } 1336 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1337 for (i = 0; i < fpdim; ++i) { 1338 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1339 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1340 npoints += Np; 1341 } 1342 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1343 for (i = 0, k = 0; i < fpdim; ++i) { 1344 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1345 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1346 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1347 } 1348 1349 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1350 PetscObject obj; 1351 PetscClassId id; 1352 PetscReal *B; 1353 PetscInt NcJ = 0, cpdim = 0, j; 1354 1355 ierr = PetscDSGetDiscretization(prob, fieldJ, &obj);CHKERRQ(ierr); 1356 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1357 if (id == PETSCFE_CLASSID) { 1358 PetscFE fe = (PetscFE) obj; 1359 1360 /* Evaluate basis at points */ 1361 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 1362 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1363 /* For now, fields only interpolate themselves */ 1364 if (fieldI == fieldJ) { 1365 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); 1366 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1367 for (i = 0, k = 0; i < fpdim; ++i) { 1368 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1369 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1370 for (p = 0; p < Np; ++p, ++k) { 1371 for (j = 0; j < cpdim; ++j) { 1372 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]; 1373 } 1374 } 1375 } 1376 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1377 } 1378 } else if (id == PETSCFV_CLASSID) { 1379 PetscFV fv = (PetscFV) obj; 1380 1381 /* Evaluate constant function at points */ 1382 ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr); 1383 cpdim = 1; 1384 /* For now, fields only interpolate themselves */ 1385 if (fieldI == fieldJ) { 1386 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); 1387 for (i = 0, k = 0; i < fpdim; ++i) { 1388 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1389 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1390 for (p = 0; p < Np; ++p, ++k) { 1391 for (j = 0; j < cpdim; ++j) { 1392 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p]; 1393 } 1394 } 1395 } 1396 } 1397 } 1398 offsetJ += cpdim*NcJ; 1399 } 1400 offsetI += fpdim*Nc; 1401 ierr = PetscFree(points);CHKERRQ(ierr); 1402 } 1403 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1404 /* Preallocate matrix */ 1405 { 1406 Mat preallocator; 1407 PetscScalar *vals; 1408 PetscInt *cellCIndices, *cellFIndices; 1409 PetscInt locRows, locCols, cell; 1410 1411 ierr = MatGetLocalSize(In, &locRows, &locCols);CHKERRQ(ierr); 1412 ierr = MatCreate(PetscObjectComm((PetscObject) In), &preallocator);CHKERRQ(ierr); 1413 ierr = MatSetType(preallocator, MATPREALLOCATOR);CHKERRQ(ierr); 1414 ierr = MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1415 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 1416 ierr = PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1417 for (cell = cStart; cell < cEnd; ++cell) { 1418 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1419 ierr = MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);CHKERRQ(ierr); 1420 } 1421 ierr = PetscFree3(vals,cellCIndices,cellFIndices);CHKERRQ(ierr); 1422 ierr = MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1423 ierr = MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1424 ierr = MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);CHKERRQ(ierr); 1425 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 1426 } 1427 /* Fill matrix */ 1428 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1429 for (c = cStart; c < cEnd; ++c) { 1430 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1431 } 1432 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1433 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1434 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1435 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1436 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1437 if (mesh->printFEM) { 1438 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1439 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1440 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1441 } 1442 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1443 PetscFunctionReturn(0); 1444 } 1445 1446 #undef __FUNCT__ 1447 #define __FUNCT__ "DMPlexComputeInterpolatorGeneral" 1448 /*@ 1449 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 1450 1451 Input Parameters: 1452 + dmf - The fine mesh 1453 . dmc - The coarse mesh 1454 - user - The user context 1455 1456 Output Parameter: 1457 . In - The interpolation matrix 1458 1459 Level: developer 1460 1461 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1462 @*/ 1463 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 1464 { 1465 DM_Plex *mesh = (DM_Plex *) dmf->data; 1466 const char *name = "Interpolator"; 1467 PetscDS prob; 1468 PetscSection fsection, csection, globalFSection, globalCSection; 1469 PetscHashJK ht; 1470 PetscLayout rLayout; 1471 PetscInt *dnz, *onz; 1472 PetscInt locRows, rStart, rEnd; 1473 PetscReal *x, *v0, *J, *invJ, detJ; 1474 PetscReal *v0c, *Jc, *invJc, detJc; 1475 PetscScalar *elemMat; 1476 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1477 PetscErrorCode ierr; 1478 1479 PetscFunctionBegin; 1480 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1481 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1482 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1483 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1484 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1485 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1486 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1487 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1488 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1489 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1490 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1491 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1492 ierr = PetscMalloc1(totDim*totDim, &elemMat);CHKERRQ(ierr); 1493 1494 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1495 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1496 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1497 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1498 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1499 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1500 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1501 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1502 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1503 for (field = 0; field < Nf; ++field) { 1504 PetscObject obj; 1505 PetscClassId id; 1506 PetscDualSpace Q = NULL; 1507 PetscQuadrature f; 1508 const PetscReal *qpoints; 1509 PetscInt Nc, Np, fpdim, i, d; 1510 1511 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1512 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1513 if (id == PETSCFE_CLASSID) { 1514 PetscFE fe = (PetscFE) obj; 1515 1516 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1517 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1518 } else if (id == PETSCFV_CLASSID) { 1519 PetscFV fv = (PetscFV) obj; 1520 1521 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1522 Nc = 1; 1523 } 1524 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1525 /* For each fine grid cell */ 1526 for (cell = cStart; cell < cEnd; ++cell) { 1527 PetscInt *findices, *cindices; 1528 PetscInt numFIndices, numCIndices; 1529 1530 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1531 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1532 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1533 for (i = 0; i < fpdim; ++i) { 1534 Vec pointVec; 1535 PetscScalar *pV; 1536 IS coarseCellIS; 1537 const PetscInt *coarseCells; 1538 PetscInt numCoarseCells, q, r, c; 1539 1540 /* Get points from the dual basis functional quadrature */ 1541 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1542 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1543 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1544 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1545 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1546 for (q = 0; q < Np; ++q) { 1547 /* Transform point to real space */ 1548 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1549 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1550 } 1551 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1552 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1553 ierr = DMLocatePoints(dmc, pointVec, &coarseCellIS);CHKERRQ(ierr); 1554 ierr = ISViewFromOptions(coarseCellIS, NULL, "-interp_is_view");CHKERRQ(ierr); 1555 /* Update preallocation info */ 1556 ierr = ISGetLocalSize(coarseCellIS, &numCoarseCells);CHKERRQ(ierr); 1557 ierr = ISGetIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1558 for (r = 0; r < Nc; ++r) { 1559 PetscHashJKKey key; 1560 PetscHashJKIter missing, iter; 1561 1562 key.j = findices[i*Nc+r]; 1563 if (key.j < 0) continue; 1564 /* Get indices for coarse elements */ 1565 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1566 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1567 for (c = 0; c < numCIndices; ++c) { 1568 key.k = cindices[c]; 1569 if (key.k < 0) continue; 1570 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1571 if (missing) { 1572 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1573 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1574 else ++onz[key.j-rStart]; 1575 } 1576 } 1577 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1578 } 1579 } 1580 ierr = ISRestoreIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1581 ierr = ISDestroy(&coarseCellIS);CHKERRQ(ierr); 1582 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1583 } 1584 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1585 } 1586 } 1587 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1588 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1589 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1590 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1591 for (field = 0; field < Nf; ++field) { 1592 PetscObject obj; 1593 PetscClassId id; 1594 PetscDualSpace Q = NULL; 1595 PetscQuadrature f; 1596 const PetscReal *qpoints, *qweights; 1597 PetscInt Nc, Np, fpdim, i, d; 1598 1599 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1600 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1601 if (id == PETSCFE_CLASSID) { 1602 PetscFE fe = (PetscFE) obj; 1603 1604 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1605 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1606 } else if (id == PETSCFV_CLASSID) { 1607 PetscFV fv = (PetscFV) obj; 1608 1609 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1610 Nc = 1; 1611 } 1612 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1613 /* For each fine grid cell */ 1614 for (cell = cStart; cell < cEnd; ++cell) { 1615 PetscInt *findices, *cindices; 1616 PetscInt numFIndices, numCIndices; 1617 1618 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1619 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1620 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1621 for (i = 0; i < fpdim; ++i) { 1622 Vec pointVec; 1623 PetscScalar *pV; 1624 IS coarseCellIS; 1625 const PetscInt *coarseCells; 1626 PetscInt numCoarseCells, cpdim, q, c, j; 1627 1628 /* Get points from the dual basis functional quadrature */ 1629 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1630 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, &qweights);CHKERRQ(ierr); 1631 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1632 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1633 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1634 for (q = 0; q < Np; ++q) { 1635 /* Transform point to real space */ 1636 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1637 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1638 } 1639 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1640 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1641 ierr = DMLocatePoints(dmc, pointVec, &coarseCellIS);CHKERRQ(ierr); 1642 /* Update preallocation info */ 1643 ierr = ISGetLocalSize(coarseCellIS, &numCoarseCells);CHKERRQ(ierr); 1644 ierr = ISGetIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1645 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1646 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1647 PetscReal pVReal[3]; 1648 1649 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1650 /* Transform points from real space to coarse reference space */ 1651 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell], NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 1652 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 1653 CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x); 1654 1655 if (id == PETSCFE_CLASSID) { 1656 PetscFE fe = (PetscFE) obj; 1657 PetscReal *B; 1658 1659 /* Evaluate coarse basis on contained point */ 1660 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1661 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 1662 /* Get elemMat entries by multiplying by weight */ 1663 for (j = 0; j < cpdim; ++j) { 1664 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = B[j*Nc + c]*qweights[ccell]; 1665 } 1666 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1667 } else { 1668 cpdim = 1; 1669 for (j = 0; j < cpdim; ++j) { 1670 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = 1.0*qweights[ccell]; 1671 } 1672 } 1673 /* Update interpolator */ 1674 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, Nc, numCIndices, elemMat);CHKERRQ(ierr);} 1675 ierr = MatSetValues(In, Nc, &findices[i*Nc], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1676 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1677 } 1678 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1679 ierr = ISRestoreIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1680 ierr = ISDestroy(&coarseCellIS);CHKERRQ(ierr); 1681 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1682 } 1683 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1684 } 1685 } 1686 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 1687 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 1688 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1689 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1690 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1691 PetscFunctionReturn(0); 1692 } 1693 1694 #undef __FUNCT__ 1695 #define __FUNCT__ "DMPlexComputeInjectorFEM" 1696 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 1697 { 1698 PetscDS prob; 1699 PetscFE *feRef; 1700 PetscFV *fvRef; 1701 Vec fv, cv; 1702 IS fis, cis; 1703 PetscSection fsection, fglobalSection, csection, cglobalSection; 1704 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 1705 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m; 1706 PetscErrorCode ierr; 1707 1708 PetscFunctionBegin; 1709 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1710 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1711 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1712 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1713 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1714 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1715 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1716 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1717 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1718 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1719 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1720 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1721 for (f = 0; f < Nf; ++f) { 1722 PetscObject obj; 1723 PetscClassId id; 1724 PetscInt fNb = 0, Nc = 0; 1725 1726 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1727 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1728 if (id == PETSCFE_CLASSID) { 1729 PetscFE fe = (PetscFE) obj; 1730 1731 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1732 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 1733 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1734 } else if (id == PETSCFV_CLASSID) { 1735 PetscFV fv = (PetscFV) obj; 1736 PetscDualSpace Q; 1737 1738 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1739 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1740 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 1741 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1742 } 1743 fTotDim += fNb*Nc; 1744 } 1745 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1746 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 1747 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 1748 PetscFE feC; 1749 PetscFV fvC; 1750 PetscDualSpace QF, QC; 1751 PetscInt NcF, NcC, fpdim, cpdim; 1752 1753 if (feRef[field]) { 1754 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 1755 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 1756 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 1757 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 1758 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1759 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 1760 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1761 } else { 1762 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 1763 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 1764 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 1765 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 1766 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1767 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 1768 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1769 } 1770 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); 1771 for (c = 0; c < cpdim; ++c) { 1772 PetscQuadrature cfunc; 1773 const PetscReal *cqpoints; 1774 PetscInt NpC; 1775 PetscBool found = PETSC_FALSE; 1776 1777 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 1778 ierr = PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 1779 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 1780 for (f = 0; f < fpdim; ++f) { 1781 PetscQuadrature ffunc; 1782 const PetscReal *fqpoints; 1783 PetscReal sum = 0.0; 1784 PetscInt NpF, comp; 1785 1786 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 1787 ierr = PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 1788 if (NpC != NpF) continue; 1789 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 1790 if (sum > 1.0e-9) continue; 1791 for (comp = 0; comp < NcC; ++comp) { 1792 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 1793 } 1794 found = PETSC_TRUE; 1795 break; 1796 } 1797 if (!found) { 1798 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 1799 if (fvRef[field]) { 1800 PetscInt comp; 1801 for (comp = 0; comp < NcC; ++comp) { 1802 cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp; 1803 } 1804 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 1805 } 1806 } 1807 offsetC += cpdim*NcC; 1808 offsetF += fpdim*NcF; 1809 } 1810 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 1811 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1812 1813 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 1814 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 1815 ierr = VecGetOwnershipRange(cv, &startC, &endC);CHKERRQ(ierr); 1816 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 1817 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 1818 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 1819 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 1820 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 1821 for (c = cStart; c < cEnd; ++c) { 1822 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 1823 for (d = 0; d < cTotDim; ++d) { 1824 if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue; 1825 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]]); 1826 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 1827 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 1828 } 1829 } 1830 ierr = PetscFree(cmap);CHKERRQ(ierr); 1831 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 1832 1833 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 1834 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 1835 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 1836 ierr = ISDestroy(&cis);CHKERRQ(ierr); 1837 ierr = ISDestroy(&fis);CHKERRQ(ierr); 1838 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 1839 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 1840 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1841 PetscFunctionReturn(0); 1842 } 1843