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