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 (f = 0; f < Nf; ++f) { 377 PetscObject obj; 378 PetscClassId id; 379 380 ierr = DMGetField(dm, f, &obj);CHKERRQ(ierr); 381 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 382 if (id == PETSCFE_CLASSID) { 383 PetscFE fe = (PetscFE) obj; 384 385 hasFE = PETSC_TRUE; 386 isFE[f] = PETSC_TRUE; 387 ierr = PetscFEGetNumComponents(fe, &numComp[f]);CHKERRQ(ierr); 388 ierr = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr); 389 } else if (id == PETSCFV_CLASSID) { 390 PetscFV fv = (PetscFV) obj; 391 392 hasFV = PETSC_TRUE; 393 isFE[f] = PETSC_FALSE; 394 ierr = PetscFVGetNumComponents(fv, &numComp[f]);CHKERRQ(ierr); 395 ierr = PetscFVGetDualSpace(fv, &cellsp[f]);CHKERRQ(ierr); 396 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 397 } 398 for (h = 0; h <= maxHeight; h++) { 399 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 400 if (!h) {pStart = cStart; pEnd = cEnd;} 401 if (pEnd <= pStart) continue; 402 totDim = 0; 403 for (f = 0; f < Nf; ++f) { 404 if (!h) { 405 sp[f] = cellsp[f]; 406 } 407 else { 408 ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);CHKERRQ(ierr); 409 if (!sp[f]) { 410 continue; 411 } 412 } 413 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 414 totDim += spDim*numComp[f]; 415 } 416 ierr = DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);CHKERRQ(ierr); 417 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim); 418 if (!totDim) continue; 419 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 420 for (p = pStart; p < pEnd; ++p) { 421 PetscFECellGeom fegeom; 422 PetscFVCellGeom fvgeom; 423 424 if (hasFE) { 425 ierr = DMPlexComputeCellGeometryFEM(dm, p, NULL, fegeom.v0, fegeom.J, NULL, &fegeom.detJ);CHKERRQ(ierr); 426 fegeom.dim = dim - h; 427 fegeom.dimEmbed = dimEmbed; 428 } 429 if (hasFV) {ierr = DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);CHKERRQ(ierr);} 430 for (f = 0, v = 0; f < Nf; ++f) { 431 void * const ctx = ctxs ? ctxs[f] : NULL; 432 433 if (!sp[f]) continue; 434 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 435 for (d = 0; d < spDim; ++d) { 436 if (funcs[f]) { 437 if (isFE[f]) {ierr = PetscDualSpaceApply(sp[f], d, time, &fegeom, numComp[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr);} 438 else {ierr = PetscDualSpaceApplyFVM(sp[f], d, time, &fvgeom, numComp[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr);} 439 } else { 440 for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0; 441 } 442 v += numComp[f]; 443 } 444 } 445 ierr = DMPlexVecSetClosure(dm, section, localX, p, values, mode);CHKERRQ(ierr); 446 } 447 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 448 if (h) { 449 for (f = 0; f < Nf; f++) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);} 450 } 451 } 452 ierr = PetscFree3(isFE, sp, numComp);CHKERRQ(ierr); 453 if (maxHeight > 0) { 454 ierr = PetscFree(cellsp);CHKERRQ(ierr); 455 } 456 PetscFunctionReturn(0); 457 } 458 459 #undef __FUNCT__ 460 #define __FUNCT__ "DMProjectFieldLocal_Plex" 461 PetscErrorCode DMProjectFieldLocal_Plex(DM dm, Vec localU, 462 void (**funcs)(PetscInt, PetscInt, PetscInt, 463 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 464 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 465 PetscReal, const PetscReal[], PetscScalar[]), 466 InsertMode mode, Vec localX) 467 { 468 DM dmAux; 469 PetscDS prob, probAux = NULL; 470 Vec A; 471 PetscSection section, sectionAux = NULL; 472 PetscDualSpace *sp; 473 PetscInt *Ncf; 474 PetscScalar *values, *u, *u_x, *a, *a_x; 475 PetscReal *x, *v0, *J, *invJ, detJ; 476 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL; 477 PetscInt Nf, NfAux = 0, dim, spDim, totDim, numValues, cStart, cEnd, cEndInterior, c, f, d, v, comp, maxHeight; 478 PetscErrorCode ierr; 479 480 PetscFunctionBegin; 481 ierr = DMPlexGetMaxProjectionHeight(dm,&maxHeight);CHKERRQ(ierr); 482 if (maxHeight > 0) {SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Field projection for height > 0 not supported yet");} 483 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 484 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 485 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 486 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 487 ierr = PetscMalloc2(Nf, &sp, Nf, &Ncf);CHKERRQ(ierr); 488 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 489 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 490 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 491 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 492 ierr = PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);CHKERRQ(ierr); 493 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 494 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 495 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 496 if (dmAux) { 497 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 498 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 499 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 500 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 501 ierr = PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);CHKERRQ(ierr); 502 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);CHKERRQ(ierr); 503 } 504 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localU, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 505 ierr = DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr); 506 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim); 507 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 508 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 509 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 510 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 511 for (c = cStart; c < cEnd; ++c) { 512 PetscScalar *coefficients = NULL, *coefficientsAux = NULL; 513 514 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 515 ierr = DMPlexVecGetClosure(dm, section, localU, c, NULL, &coefficients);CHKERRQ(ierr); 516 if (dmAux) {ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);CHKERRQ(ierr);} 517 for (f = 0, v = 0; f < Nf; ++f) { 518 PetscObject obj; 519 PetscClassId id; 520 521 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 522 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 523 if (id == PETSCFE_CLASSID) { 524 PetscFE fe = (PetscFE) obj; 525 526 ierr = PetscFEGetDualSpace(fe, &sp[f]);CHKERRQ(ierr); 527 ierr = PetscFEGetNumComponents(fe, &Ncf[f]);CHKERRQ(ierr); 528 } else if (id == PETSCFV_CLASSID) { 529 PetscFV fv = (PetscFV) obj; 530 531 ierr = PetscFVGetNumComponents(fv, &Ncf[f]);CHKERRQ(ierr); 532 ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr); 533 } 534 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 535 for (d = 0; d < spDim; ++d) { 536 PetscQuadrature quad; 537 const PetscReal *points, *weights; 538 PetscInt numPoints, q; 539 540 if (funcs[f]) { 541 ierr = PetscDualSpaceGetFunctional(sp[f], d, &quad);CHKERRQ(ierr); 542 ierr = PetscQuadratureGetData(quad, NULL, &numPoints, &points, &weights);CHKERRQ(ierr); 543 for (q = 0; q < numPoints; ++q) { 544 CoordinatesRefToReal(dim, dim, v0, J, &points[q*dim], x); 545 ierr = EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, coefficients, NULL, u, u_x, NULL);CHKERRQ(ierr); 546 ierr = EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, coefficientsAux, NULL, a, a_x, NULL);CHKERRQ(ierr); 547 (*funcs[f])(dim, Nf, NfAux, uOff, uOff_x, u, NULL, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, x, &values[v]); 548 } 549 } else { 550 for (comp = 0; comp < Ncf[f]; ++comp) values[v+comp] = 0.0; 551 } 552 v += Ncf[f]; 553 } 554 } 555 ierr = DMPlexVecRestoreClosure(dm, section, localU, c, NULL, &coefficients);CHKERRQ(ierr); 556 if (dmAux) {ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);CHKERRQ(ierr);} 557 ierr = DMPlexVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr); 558 } 559 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 560 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 561 ierr = PetscFree2(sp, Ncf);CHKERRQ(ierr); 562 PetscFunctionReturn(0); 563 } 564 565 #undef __FUNCT__ 566 #define __FUNCT__ "DMPlexInsertBoundaryValues_FEM_Internal" 567 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) 568 { 569 PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx); 570 void **ctxs; 571 PetscInt numFields; 572 PetscErrorCode ierr; 573 574 PetscFunctionBegin; 575 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 576 ierr = PetscCalloc2(numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 577 funcs[field] = func; 578 ctxs[field] = ctx; 579 ierr = DMProjectFunctionLabelLocal(dm, time, label, numids, ids, funcs, ctxs, INSERT_BC_VALUES, locX);CHKERRQ(ierr); 580 ierr = PetscFree2(funcs,ctxs);CHKERRQ(ierr); 581 PetscFunctionReturn(0); 582 } 583 584 #undef __FUNCT__ 585 #define __FUNCT__ "DMPlexInsertBoundaryValues_FVM_Internal" 586 /* This ignores numcomps/comps */ 587 static PetscErrorCode DMPlexInsertBoundaryValues_FVM_Internal(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad, 588 PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX) 589 { 590 PetscDS prob; 591 PetscSF sf; 592 DM dmFace, dmCell, dmGrad; 593 const PetscScalar *facegeom, *cellgeom = NULL, *grad; 594 const PetscInt *leaves; 595 PetscScalar *x, *fx; 596 PetscInt dim, nleaves, loc, fStart, fEnd, pdim, i; 597 PetscErrorCode ierr; 598 599 PetscFunctionBegin; 600 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 601 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 602 nleaves = PetscMax(0, nleaves); 603 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 604 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 605 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 606 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 607 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 608 if (cellGeometry) { 609 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 610 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 611 } 612 if (Grad) { 613 PetscFV fv; 614 615 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);CHKERRQ(ierr); 616 ierr = VecGetDM(Grad, &dmGrad);CHKERRQ(ierr); 617 ierr = VecGetArrayRead(Grad, &grad);CHKERRQ(ierr); 618 ierr = PetscFVGetNumComponents(fv, &pdim);CHKERRQ(ierr); 619 ierr = DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 620 } 621 ierr = VecGetArray(locX, &x);CHKERRQ(ierr); 622 for (i = 0; i < numids; ++i) { 623 IS faceIS; 624 const PetscInt *faces; 625 PetscInt numFaces, f; 626 627 ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr); 628 if (!faceIS) continue; /* No points with that id on this process */ 629 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 630 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 631 for (f = 0; f < numFaces; ++f) { 632 const PetscInt face = faces[f], *cells; 633 PetscFVFaceGeom *fg; 634 635 if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */ 636 ierr = PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);CHKERRQ(ierr); 637 if (loc >= 0) continue; 638 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 639 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 640 if (Grad) { 641 PetscFVCellGeom *cg; 642 PetscScalar *cx, *cgrad; 643 PetscScalar *xG; 644 PetscReal dx[3]; 645 PetscInt d; 646 647 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);CHKERRQ(ierr); 648 ierr = DMPlexPointLocalRead(dm, cells[0], x, &cx);CHKERRQ(ierr); 649 ierr = DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);CHKERRQ(ierr); 650 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 651 DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx); 652 for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx); 653 ierr = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);CHKERRQ(ierr); 654 } else { 655 PetscScalar *xI; 656 PetscScalar *xG; 657 658 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr); 659 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 660 ierr = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);CHKERRQ(ierr); 661 } 662 } 663 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 664 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 665 } 666 ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr); 667 if (Grad) { 668 ierr = DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 669 ierr = VecRestoreArrayRead(Grad, &grad);CHKERRQ(ierr); 670 } 671 if (cellGeometry) { 672 ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 673 } 674 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 675 PetscFunctionReturn(0); 676 } 677 678 #undef __FUNCT__ 679 #define __FUNCT__ "DMPlexInsertBoundaryValues" 680 PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM) 681 { 682 PetscInt numBd, b; 683 PetscErrorCode ierr; 684 685 PetscFunctionBegin; 686 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 687 PetscValidHeaderSpecific(locX, VEC_CLASSID, 2); 688 if (faceGeomFVM) {PetscValidHeaderSpecific(faceGeomFVM, VEC_CLASSID, 4);} 689 if (cellGeomFVM) {PetscValidHeaderSpecific(cellGeomFVM, VEC_CLASSID, 5);} 690 if (gradFVM) {PetscValidHeaderSpecific(gradFVM, VEC_CLASSID, 6);} 691 ierr = DMGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 692 for (b = 0; b < numBd; ++b) { 693 PetscBool isEssential; 694 const char *labelname; 695 DMLabel label; 696 PetscInt field; 697 PetscObject obj; 698 PetscClassId id; 699 void (*func)(); 700 PetscInt numids; 701 const PetscInt *ids; 702 void *ctx; 703 704 ierr = DMGetBoundary(dm, b, &isEssential, NULL, &labelname, &field, NULL, NULL, &func, &numids, &ids, &ctx);CHKERRQ(ierr); 705 if (insertEssential != isEssential) continue; 706 ierr = DMGetLabel(dm, labelname, &label);CHKERRQ(ierr); 707 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 708 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 709 if (id == PETSCFE_CLASSID) { 710 if (!isEssential) continue; /* for FEM, there is no insertion to be done for non-essential boundary conditions */ 711 ierr = DMPlexLabelAddCells(dm,label);CHKERRQ(ierr); 712 ierr = DMPlexInsertBoundaryValues_FEM_Internal(dm, time, field, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);CHKERRQ(ierr); 713 ierr = DMPlexLabelClearCells(dm,label);CHKERRQ(ierr); 714 } else if (id == PETSCFV_CLASSID) { 715 if (!faceGeomFVM) continue; 716 ierr = DMPlexInsertBoundaryValues_FVM_Internal(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, 717 field, label, numids, ids, (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);CHKERRQ(ierr); 718 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 719 } 720 PetscFunctionReturn(0); 721 } 722 723 #undef __FUNCT__ 724 #define __FUNCT__ "DMComputeL2Diff_Plex" 725 PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 726 { 727 const PetscInt debug = 0; 728 PetscSection section; 729 PetscQuadrature quad; 730 Vec localX; 731 PetscScalar *funcVal, *interpolant; 732 PetscReal *coords, *v0, *J, *invJ, detJ; 733 PetscReal localDiff = 0.0; 734 const PetscReal *quadPoints, *quadWeights; 735 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 736 PetscErrorCode ierr; 737 738 PetscFunctionBegin; 739 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 740 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 741 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 742 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 743 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 744 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 745 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 746 for (field = 0; field < numFields; ++field) { 747 PetscObject obj; 748 PetscClassId id; 749 PetscInt Nc; 750 751 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 752 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 753 if (id == PETSCFE_CLASSID) { 754 PetscFE fe = (PetscFE) obj; 755 756 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 757 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 758 } else if (id == PETSCFV_CLASSID) { 759 PetscFV fv = (PetscFV) obj; 760 761 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 762 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 763 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 764 numComponents += Nc; 765 } 766 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 767 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 768 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 769 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 770 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 771 for (c = cStart; c < cEnd; ++c) { 772 PetscScalar *x = NULL; 773 PetscReal elemDiff = 0.0; 774 775 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 776 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 777 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 778 779 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 780 PetscObject obj; 781 PetscClassId id; 782 void * const ctx = ctxs ? ctxs[field] : NULL; 783 PetscInt Nb, Nc, q, fc; 784 785 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 786 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 787 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 788 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 789 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 790 if (debug) { 791 char title[1024]; 792 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 793 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 794 } 795 for (q = 0; q < numQuadPoints; ++q) { 796 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 797 ierr = (*funcs[field])(dim, time, coords, Nc, funcVal, ctx);CHKERRQ(ierr); 798 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 799 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 800 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 801 for (fc = 0; fc < Nc; ++fc) { 802 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);} 803 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 804 } 805 } 806 fieldOffset += Nb*Nc; 807 } 808 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 809 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 810 localDiff += elemDiff; 811 } 812 ierr = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 813 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 814 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 815 *diff = PetscSqrtReal(*diff); 816 PetscFunctionReturn(0); 817 } 818 819 #undef __FUNCT__ 820 #define __FUNCT__ "DMComputeL2GradientDiff_Plex" 821 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) 822 { 823 const PetscInt debug = 0; 824 PetscSection section; 825 PetscQuadrature quad; 826 Vec localX; 827 PetscScalar *funcVal, *interpolantVec; 828 PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ; 829 PetscReal localDiff = 0.0; 830 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset, comp; 831 PetscErrorCode ierr; 832 833 PetscFunctionBegin; 834 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 835 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 836 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 837 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 838 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 839 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 840 for (field = 0; field < numFields; ++field) { 841 PetscFE fe; 842 PetscInt Nc; 843 844 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 845 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 846 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 847 numComponents += Nc; 848 } 849 /* ierr = DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 850 ierr = PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr); 851 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 852 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 853 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 854 for (c = cStart; c < cEnd; ++c) { 855 PetscScalar *x = NULL; 856 PetscReal elemDiff = 0.0; 857 858 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 859 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 860 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 861 862 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 863 PetscFE fe; 864 void * const ctx = ctxs ? ctxs[field] : NULL; 865 const PetscReal *quadPoints, *quadWeights; 866 PetscReal *basisDer; 867 PetscInt numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g; 868 869 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 870 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 871 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 872 ierr = PetscFEGetNumComponents(fe, &Ncomp);CHKERRQ(ierr); 873 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 874 if (debug) { 875 char title[1024]; 876 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 877 ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr); 878 } 879 for (q = 0; q < numQuadPoints; ++q) { 880 for (d = 0; d < dim; d++) { 881 coords[d] = v0[d]; 882 for (e = 0; e < dim; e++) { 883 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 884 } 885 } 886 ierr = (*funcs[field])(dim, time, coords, n, numFields, funcVal, ctx);CHKERRQ(ierr); 887 for (fc = 0; fc < Ncomp; ++fc) { 888 PetscScalar interpolant = 0.0; 889 890 for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0; 891 for (f = 0; f < Nb; ++f) { 892 const PetscInt fidx = f*Ncomp+fc; 893 894 for (d = 0; d < dim; ++d) { 895 realSpaceDer[d] = 0.0; 896 for (g = 0; g < dim; ++g) { 897 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g]; 898 } 899 interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d]; 900 } 901 } 902 for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d]; 903 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);} 904 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 905 } 906 } 907 comp += Ncomp; 908 fieldOffset += Nb*Ncomp; 909 } 910 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 911 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 912 localDiff += elemDiff; 913 } 914 ierr = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr); 915 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 916 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 917 *diff = PetscSqrtReal(*diff); 918 PetscFunctionReturn(0); 919 } 920 921 #undef __FUNCT__ 922 #define __FUNCT__ "DMComputeL2FieldDiff_Plex" 923 PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 924 { 925 const PetscInt debug = 0; 926 PetscSection section; 927 PetscQuadrature quad; 928 Vec localX; 929 PetscScalar *funcVal, *interpolant; 930 PetscReal *coords, *v0, *J, *invJ, detJ; 931 PetscReal *localDiff; 932 const PetscReal *quadPoints, *quadWeights; 933 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 934 PetscErrorCode ierr; 935 936 PetscFunctionBegin; 937 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 938 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 939 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 940 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 941 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 942 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 943 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 944 for (field = 0; field < numFields; ++field) { 945 PetscObject obj; 946 PetscClassId id; 947 PetscInt Nc; 948 949 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 950 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 951 if (id == PETSCFE_CLASSID) { 952 PetscFE fe = (PetscFE) obj; 953 954 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 955 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 956 } else if (id == PETSCFV_CLASSID) { 957 PetscFV fv = (PetscFV) obj; 958 959 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 960 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 961 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 962 numComponents += Nc; 963 } 964 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 965 ierr = PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 966 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 967 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 968 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 969 for (c = cStart; c < cEnd; ++c) { 970 PetscScalar *x = NULL; 971 972 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 973 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 974 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 975 976 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 977 PetscObject obj; 978 PetscClassId id; 979 void * const ctx = ctxs ? ctxs[field] : NULL; 980 PetscInt Nb, Nc, q, fc; 981 982 PetscReal elemDiff = 0.0; 983 984 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 985 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 986 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 987 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 988 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 989 if (debug) { 990 char title[1024]; 991 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 992 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 993 } 994 for (q = 0; q < numQuadPoints; ++q) { 995 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 996 ierr = (*funcs[field])(dim, time, coords, numFields, funcVal, ctx);CHKERRQ(ierr); 997 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 998 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 999 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1000 for (fc = 0; fc < Nc; ++fc) { 1001 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);} 1002 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 1003 } 1004 } 1005 fieldOffset += Nb*Nc; 1006 localDiff[field] += elemDiff; 1007 } 1008 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1009 } 1010 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1011 ierr = MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1012 for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]); 1013 ierr = PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 1014 PetscFunctionReturn(0); 1015 } 1016 1017 #undef __FUNCT__ 1018 #define __FUNCT__ "DMPlexComputeL2DiffVec" 1019 /*@C 1020 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. 1021 1022 Input Parameters: 1023 + dm - The DM 1024 . time - The time 1025 . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation 1026 . ctxs - Optional array of contexts to pass to each function, or NULL. 1027 - X - The coefficient vector u_h 1028 1029 Output Parameter: 1030 . D - A Vec which holds the difference ||u - u_h||_2 for each cell 1031 1032 Level: developer 1033 1034 .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff() 1035 @*/ 1036 PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D) 1037 { 1038 PetscSection section; 1039 PetscQuadrature quad; 1040 Vec localX; 1041 PetscScalar *funcVal, *interpolant; 1042 PetscReal *coords, *v0, *J, *invJ, detJ; 1043 const PetscReal *quadPoints, *quadWeights; 1044 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 1045 PetscErrorCode ierr; 1046 1047 PetscFunctionBegin; 1048 ierr = VecSet(D, 0.0);CHKERRQ(ierr); 1049 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1050 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1051 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1052 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1053 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 1054 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1055 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1056 for (field = 0; field < numFields; ++field) { 1057 PetscObject obj; 1058 PetscClassId id; 1059 PetscInt Nc; 1060 1061 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1062 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1063 if (id == PETSCFE_CLASSID) { 1064 PetscFE fe = (PetscFE) obj; 1065 1066 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1067 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1068 } else if (id == PETSCFV_CLASSID) { 1069 PetscFV fv = (PetscFV) obj; 1070 1071 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1072 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1073 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1074 numComponents += Nc; 1075 } 1076 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 1077 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1078 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1079 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1080 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1081 for (c = cStart; c < cEnd; ++c) { 1082 PetscScalar *x = NULL; 1083 PetscScalar elemDiff = 0.0; 1084 1085 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1086 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 1087 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1088 1089 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1090 PetscObject obj; 1091 PetscClassId id; 1092 void * const ctx = ctxs ? ctxs[field] : NULL; 1093 PetscInt Nb, Nc, q, fc; 1094 1095 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1096 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1097 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1098 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1099 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1100 if (funcs[field]) { 1101 for (q = 0; q < numQuadPoints; ++q) { 1102 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 1103 ierr = (*funcs[field])(dim, time, coords, Nc, funcVal, ctx);CHKERRQ(ierr); 1104 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1105 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1106 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1107 for (fc = 0; fc < Nc; ++fc) { 1108 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 1109 } 1110 } 1111 } 1112 fieldOffset += Nb*Nc; 1113 } 1114 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1115 ierr = VecSetValue(D, c - cStart, elemDiff, INSERT_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 = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 1158 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1159 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1160 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1161 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1162 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1163 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1164 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1165 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1166 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1167 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1168 numCells = cEnd - cStart; 1169 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1170 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1171 if (dmAux) { 1172 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1173 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1174 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1175 } 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, NULL);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, NULL);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, NULL);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, NULL);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, NULL);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, NULL);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, NULL);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, NULL);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