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