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