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