1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petscsf.h> 3 4 #include <petsc/private/petscfeimpl.h> 5 #include <petsc/private/petscfvimpl.h> 6 7 #undef __FUNCT__ 8 #define __FUNCT__ "DMPlexGetScale" 9 PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale) 10 { 11 DM_Plex *mesh = (DM_Plex*) dm->data; 12 13 PetscFunctionBegin; 14 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 15 PetscValidPointer(scale, 3); 16 *scale = mesh->scale[unit]; 17 PetscFunctionReturn(0); 18 } 19 20 #undef __FUNCT__ 21 #define __FUNCT__ "DMPlexSetScale" 22 PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale) 23 { 24 DM_Plex *mesh = (DM_Plex*) dm->data; 25 26 PetscFunctionBegin; 27 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 28 mesh->scale[unit] = scale; 29 PetscFunctionReturn(0); 30 } 31 32 PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k) 33 { 34 switch (i) { 35 case 0: 36 switch (j) { 37 case 0: return 0; 38 case 1: 39 switch (k) { 40 case 0: return 0; 41 case 1: return 0; 42 case 2: return 1; 43 } 44 case 2: 45 switch (k) { 46 case 0: return 0; 47 case 1: return -1; 48 case 2: return 0; 49 } 50 } 51 case 1: 52 switch (j) { 53 case 0: 54 switch (k) { 55 case 0: return 0; 56 case 1: return 0; 57 case 2: return -1; 58 } 59 case 1: return 0; 60 case 2: 61 switch (k) { 62 case 0: return 1; 63 case 1: return 0; 64 case 2: return 0; 65 } 66 } 67 case 2: 68 switch (j) { 69 case 0: 70 switch (k) { 71 case 0: return 0; 72 case 1: return 1; 73 case 2: return 0; 74 } 75 case 1: 76 switch (k) { 77 case 0: return -1; 78 case 1: return 0; 79 case 2: return 0; 80 } 81 case 2: return 0; 82 } 83 } 84 return 0; 85 } 86 87 #undef __FUNCT__ 88 #define __FUNCT__ "DMPlexProjectRigidBody" 89 static PetscErrorCode DMPlexProjectRigidBody(PetscInt dim, PetscReal t, const PetscReal X[], PetscInt Nf, PetscScalar *mode, void *ctx) 90 { 91 PetscInt *ctxInt = (PetscInt *) ctx; 92 PetscInt dim2 = ctxInt[0]; 93 PetscInt d = ctxInt[1]; 94 PetscInt i, j, k = dim > 2 ? d - dim : d; 95 96 PetscFunctionBegin; 97 if (dim != dim2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input dimension %d does not match context dimension %d", dim, dim2); 98 for (i = 0; i < dim; i++) mode[i] = 0.; 99 if (d < dim) { 100 mode[d] = 1.; 101 } else { 102 for (i = 0; i < dim; i++) { 103 for (j = 0; j < dim; j++) { 104 mode[j] += epsilon(i, j, k)*X[i]; 105 } 106 } 107 } 108 PetscFunctionReturn(0); 109 } 110 111 #undef __FUNCT__ 112 #define __FUNCT__ "DMPlexCreateRigidBody" 113 /*@C 114 DMPlexCreateRigidBody - for the default global section, create rigid body modes from coordinates 115 116 Collective on DM 117 118 Input Arguments: 119 . dm - the DM 120 121 Output Argument: 122 . sp - the null space 123 124 Note: This is necessary to take account of Dirichlet conditions on the displacements 125 126 Level: advanced 127 128 .seealso: MatNullSpaceCreate() 129 @*/ 130 PetscErrorCode DMPlexCreateRigidBody(DM dm, MatNullSpace *sp) 131 { 132 MPI_Comm comm; 133 Vec mode[6]; 134 PetscSection section, globalSection; 135 PetscInt dim, dimEmbed, n, m, d, i, j; 136 PetscErrorCode ierr; 137 138 PetscFunctionBegin; 139 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 140 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 141 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 142 if (dim == 1) { 143 ierr = MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);CHKERRQ(ierr); 144 PetscFunctionReturn(0); 145 } 146 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 147 ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr); 148 ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr); 149 m = (dim*(dim+1))/2; 150 ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr); 151 ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr); 152 ierr = VecSetUp(mode[0]);CHKERRQ(ierr); 153 for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);} 154 for (d = 0; d < m; d++) { 155 PetscInt ctx[2]; 156 void *voidctx = (void *) (&ctx[0]); 157 PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody; 158 159 ctx[0] = dimEmbed; 160 ctx[1] = d; 161 ierr = DMProjectFunction(dm, 0.0, &func, &voidctx, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 162 } 163 for (i = 0; i < dim; ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);} 164 /* Orthonormalize system */ 165 for (i = dim; i < m; ++i) { 166 PetscScalar dots[6]; 167 168 ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr); 169 for (j = 0; j < i; ++j) dots[j] *= -1.0; 170 ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr); 171 ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr); 172 } 173 ierr = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);CHKERRQ(ierr); 174 for (i = 0; i< m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);} 175 PetscFunctionReturn(0); 176 } 177 178 #undef __FUNCT__ 179 #define __FUNCT__ "DMPlexSetMaxProjectionHeight" 180 /*@ 181 DMPlexSetMaxProjectionHeight - In DMPlexProjectXXXLocal() functions, the projected values of a basis function's dofs 182 are computed by associating the basis function with one of the mesh points in its transitively-closed support, and 183 evaluating the dual space basis of that point. A basis function is associated with the point in its 184 transitively-closed support whose mesh height is highest (w.r.t. DAG height), but not greater than the maximum 185 projection height, which is set with this function. By default, the maximum projection height is zero, which means 186 that only mesh cells are used to project basis functions. A height of one, for example, evaluates a cell-interior 187 basis functions using its cells dual space basis, but all other basis functions with the dual space basis of a face. 188 189 Input Parameters: 190 + dm - the DMPlex object 191 - height - the maximum projection height >= 0 192 193 Level: advanced 194 195 .seealso: DMPlexGetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal() 196 @*/ 197 PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height) 198 { 199 DM_Plex *plex = (DM_Plex *) dm->data; 200 201 PetscFunctionBegin; 202 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 203 plex->maxProjectionHeight = height; 204 PetscFunctionReturn(0); 205 } 206 207 #undef __FUNCT__ 208 #define __FUNCT__ "DMPlexGetMaxProjectionHeight" 209 /*@ 210 DMPlexGetMaxProjectionHeight - Get the maximum height (w.r.t. DAG) of mesh points used to evaluate dual bases in 211 DMPlexProjectXXXLocal() functions. 212 213 Input Parameters: 214 . dm - the DMPlex object 215 216 Output Parameters: 217 . height - the maximum projection height 218 219 Level: intermediate 220 221 .seealso: DMPlexSetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal() 222 @*/ 223 PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height) 224 { 225 DM_Plex *plex = (DM_Plex *) dm->data; 226 227 PetscFunctionBegin; 228 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 229 *height = plex->maxProjectionHeight; 230 PetscFunctionReturn(0); 231 } 232 233 #undef __FUNCT__ 234 #define __FUNCT__ "DMProjectFunctionLabelLocal_Plex" 235 PetscErrorCode DMProjectFunctionLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 236 { 237 PetscDualSpace *sp, *cellsp; 238 PetscInt *numComp; 239 PetscSection section; 240 PetscScalar *values; 241 PetscBool *fieldActive; 242 PetscInt numFields, dim, dimEmbed, spDim, totDim = 0, numValues, pStart, pEnd, cStart, cEnd, cEndInterior, f, d, v, i, comp, maxHeight, h; 243 PetscErrorCode ierr; 244 245 PetscFunctionBegin; 246 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 247 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 248 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 249 if (cEnd <= cStart) PetscFunctionReturn(0); 250 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 251 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 252 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 253 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 254 ierr = PetscMalloc2(numFields,&sp,numFields,&numComp);CHKERRQ(ierr); 255 ierr = DMPlexGetMaxProjectionHeight(dm,&maxHeight);CHKERRQ(ierr); 256 if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"maximum projection height %d not in [0, %d)\n", maxHeight,dim);} 257 if (maxHeight > 0) {ierr = PetscMalloc1(numFields,&cellsp);CHKERRQ(ierr);} 258 else {cellsp = sp;} 259 for (h = 0; h <= maxHeight; h++) { 260 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 261 if (!h) {pStart = cStart; pEnd = cEnd;} 262 if (pEnd <= pStart) continue; 263 totDim = 0; 264 for (f = 0; f < numFields; ++f) { 265 PetscObject obj; 266 PetscClassId id; 267 268 ierr = DMGetField(dm, f, &obj);CHKERRQ(ierr); 269 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 270 if (id == PETSCFE_CLASSID) { 271 PetscFE fe = (PetscFE) obj; 272 273 ierr = PetscFEGetNumComponents(fe, &numComp[f]);CHKERRQ(ierr); 274 if (!h) { 275 ierr = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr); 276 sp[f] = cellsp[f]; 277 ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr); 278 } else { 279 ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);CHKERRQ(ierr); 280 if (!sp[f]) continue; 281 } 282 } else if (id == PETSCFV_CLASSID) { 283 PetscFV fv = (PetscFV) obj; 284 285 ierr = PetscFVGetNumComponents(fv, &numComp[f]);CHKERRQ(ierr); 286 ierr = PetscFVGetDualSpace(fv, &sp[f]);CHKERRQ(ierr); 287 ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr); 288 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 289 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 290 totDim += spDim*numComp[f]; 291 } 292 ierr = DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);CHKERRQ(ierr); 293 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim); 294 if (!totDim) continue; 295 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 296 ierr = DMGetWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);CHKERRQ(ierr); 297 for (f = 0; f < numFields; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE; 298 for (i = 0; i < numIds; ++i) { 299 IS pointIS; 300 const PetscInt *points; 301 PetscInt n, p; 302 303 ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr); 304 ierr = ISGetLocalSize(pointIS, &n);CHKERRQ(ierr); 305 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 306 for (p = 0; p < n; ++p) { 307 const PetscInt point = points[p]; 308 PetscFECellGeom geom; 309 310 if ((point < pStart) || (point >= pEnd)) continue; 311 ierr = DMPlexComputeCellGeometryFEM(dm, point, NULL, geom.v0, geom.J, NULL, &geom.detJ);CHKERRQ(ierr); 312 geom.dim = dim - h; 313 geom.dimEmbed = dimEmbed; 314 for (f = 0, v = 0; f < numFields; ++f) { 315 void * const ctx = ctxs ? ctxs[f] : NULL; 316 317 if (!sp[f]) continue; 318 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 319 for (d = 0; d < spDim; ++d) { 320 if (funcs[f]) { 321 ierr = PetscDualSpaceApply(sp[f], d, time, &geom, numComp[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr); 322 } else { 323 for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0; 324 } 325 v += numComp[f]; 326 } 327 } 328 ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, values, mode);CHKERRQ(ierr); 329 } 330 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 331 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 332 } 333 if (h) { 334 for (f = 0; f < numFields; ++f) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);} 335 } 336 } 337 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 338 ierr = DMRestoreWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);CHKERRQ(ierr); 339 for (f = 0; f < numFields; ++f) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);} 340 ierr = PetscFree2(sp, numComp);CHKERRQ(ierr); 341 if (maxHeight > 0) { 342 ierr = PetscFree(cellsp);CHKERRQ(ierr); 343 } 344 PetscFunctionReturn(0); 345 } 346 347 #undef __FUNCT__ 348 #define __FUNCT__ "DMProjectFunctionLocal_Plex" 349 PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 350 { 351 PetscDualSpace *sp, *cellsp; 352 PetscInt *numComp; 353 PetscSection section; 354 PetscScalar *values; 355 PetscInt Nf, dim, dimEmbed, spDim, totDim = 0, numValues, pStart, pEnd, p, cStart, cEnd, cEndInterior, f, d, v, comp, h, maxHeight; 356 PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE; 357 PetscErrorCode ierr; 358 359 PetscFunctionBegin; 360 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 361 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 362 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 363 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 364 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 365 ierr = PetscMalloc3(Nf, &isFE, Nf, &sp, Nf, &numComp);CHKERRQ(ierr); 366 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 367 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 368 ierr = DMPlexGetMaxProjectionHeight(dm,&maxHeight);CHKERRQ(ierr); 369 if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"maximum projection height %d not in [0, %d)\n", maxHeight,dim);} 370 if (maxHeight > 0) { 371 ierr = PetscMalloc1(Nf,&cellsp);CHKERRQ(ierr); 372 } 373 else { 374 cellsp = sp; 375 } 376 for (f = 0; f < Nf; ++f) { 377 PetscObject obj; 378 PetscClassId id; 379 380 ierr = DMGetField(dm, f, &obj);CHKERRQ(ierr); 381 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 382 if (id == PETSCFE_CLASSID) { 383 PetscFE fe = (PetscFE) obj; 384 385 hasFE = PETSC_TRUE; 386 isFE[f] = PETSC_TRUE; 387 ierr = PetscFEGetNumComponents(fe, &numComp[f]);CHKERRQ(ierr); 388 ierr = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr); 389 } else if (id == PETSCFV_CLASSID) { 390 PetscFV fv = (PetscFV) obj; 391 392 hasFV = PETSC_TRUE; 393 isFE[f] = PETSC_FALSE; 394 ierr = PetscFVGetNumComponents(fv, &numComp[f]);CHKERRQ(ierr); 395 ierr = PetscFVGetDualSpace(fv, &cellsp[f]);CHKERRQ(ierr); 396 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 397 } 398 for (h = 0; h <= maxHeight; h++) { 399 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 400 if (!h) {pStart = cStart; pEnd = cEnd;} 401 if (pEnd <= pStart) continue; 402 totDim = 0; 403 for (f = 0; f < Nf; ++f) { 404 if (!h) { 405 sp[f] = cellsp[f]; 406 } 407 else { 408 ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);CHKERRQ(ierr); 409 if (!sp[f]) { 410 continue; 411 } 412 } 413 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 414 totDim += spDim*numComp[f]; 415 } 416 ierr = DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);CHKERRQ(ierr); 417 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim); 418 if (!totDim) continue; 419 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 420 for (p = pStart; p < pEnd; ++p) { 421 PetscFECellGeom fegeom; 422 PetscFVCellGeom fvgeom; 423 424 if (hasFE) { 425 ierr = DMPlexComputeCellGeometryFEM(dm, p, NULL, fegeom.v0, fegeom.J, NULL, &fegeom.detJ);CHKERRQ(ierr); 426 fegeom.dim = dim - h; 427 fegeom.dimEmbed = dimEmbed; 428 } 429 if (hasFV) {ierr = DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);CHKERRQ(ierr);} 430 for (f = 0, v = 0; f < Nf; ++f) { 431 void * const ctx = ctxs ? ctxs[f] : NULL; 432 433 if (!sp[f]) continue; 434 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 435 for (d = 0; d < spDim; ++d) { 436 if (funcs[f]) { 437 if (isFE[f]) {ierr = PetscDualSpaceApply(sp[f], d, time, &fegeom, numComp[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr);} 438 else {ierr = PetscDualSpaceApplyFVM(sp[f], d, time, &fvgeom, numComp[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr);} 439 } else { 440 for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0; 441 } 442 v += numComp[f]; 443 } 444 } 445 ierr = DMPlexVecSetClosure(dm, section, localX, p, values, mode);CHKERRQ(ierr); 446 } 447 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 448 if (h) { 449 for (f = 0; f < Nf; f++) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);} 450 } 451 } 452 ierr = PetscFree3(isFE, sp, numComp);CHKERRQ(ierr); 453 if (maxHeight > 0) { 454 ierr = PetscFree(cellsp);CHKERRQ(ierr); 455 } 456 PetscFunctionReturn(0); 457 } 458 459 #undef __FUNCT__ 460 #define __FUNCT__ "DMProjectFieldLocal_Plex" 461 PetscErrorCode DMProjectFieldLocal_Plex(DM dm, Vec localU, 462 void (**funcs)(PetscInt, PetscInt, PetscInt, 463 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 464 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 465 PetscReal, const PetscReal[], PetscScalar[]), 466 InsertMode mode, Vec localX) 467 { 468 DM dmAux; 469 PetscDS prob, probAux = NULL; 470 Vec A; 471 PetscSection section, sectionAux = NULL; 472 PetscDualSpace *sp; 473 PetscInt *Ncf; 474 PetscScalar *values, *u, *u_x, *a, *a_x; 475 PetscReal *x, *v0, *J, *invJ, detJ; 476 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL; 477 PetscInt Nf, NfAux = 0, dim, spDim, totDim, numValues, cStart, cEnd, cEndInterior, c, f, d, v, comp, maxHeight; 478 PetscErrorCode ierr; 479 480 PetscFunctionBegin; 481 ierr = DMPlexGetMaxProjectionHeight(dm,&maxHeight);CHKERRQ(ierr); 482 if (maxHeight > 0) {SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Field projection for height > 0 not supported yet");} 483 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 484 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 485 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 486 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 487 ierr = PetscMalloc2(Nf, &sp, Nf, &Ncf);CHKERRQ(ierr); 488 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 489 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 490 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 491 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 492 ierr = PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);CHKERRQ(ierr); 493 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 494 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 495 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 496 if (dmAux) { 497 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 498 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 499 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 500 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 501 ierr = PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);CHKERRQ(ierr); 502 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);CHKERRQ(ierr); 503 } 504 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localU, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 505 ierr = DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr); 506 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim); 507 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 508 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 509 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 510 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 511 for (c = cStart; c < cEnd; ++c) { 512 PetscScalar *coefficients = NULL, *coefficientsAux = NULL; 513 514 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 515 ierr = DMPlexVecGetClosure(dm, section, localU, c, NULL, &coefficients);CHKERRQ(ierr); 516 if (dmAux) {ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);CHKERRQ(ierr);} 517 for (f = 0, v = 0; f < Nf; ++f) { 518 PetscObject obj; 519 PetscClassId id; 520 521 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 522 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 523 if (id == PETSCFE_CLASSID) { 524 PetscFE fe = (PetscFE) obj; 525 526 ierr = PetscFEGetDualSpace(fe, &sp[f]);CHKERRQ(ierr); 527 ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr); 528 ierr = PetscFEGetNumComponents(fe, &Ncf[f]);CHKERRQ(ierr); 529 } else if (id == PETSCFV_CLASSID) { 530 PetscFV fv = (PetscFV) obj; 531 532 ierr = PetscFVGetNumComponents(fv, &Ncf[f]);CHKERRQ(ierr); 533 ierr = PetscFVGetDualSpace(fv, &sp[f]);CHKERRQ(ierr); 534 ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr); 535 } 536 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 537 for (d = 0; d < spDim; ++d) { 538 PetscQuadrature quad; 539 const PetscReal *points, *weights; 540 PetscInt numPoints, q; 541 542 if (funcs[f]) { 543 ierr = PetscDualSpaceGetFunctional(sp[f], d, &quad);CHKERRQ(ierr); 544 ierr = PetscQuadratureGetData(quad, NULL, &numPoints, &points, &weights);CHKERRQ(ierr); 545 for (q = 0; q < numPoints; ++q) { 546 CoordinatesRefToReal(dim, dim, v0, J, &points[q*dim], x); 547 ierr = EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, coefficients, NULL, u, u_x, NULL);CHKERRQ(ierr); 548 ierr = EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, coefficientsAux, NULL, a, a_x, NULL);CHKERRQ(ierr); 549 (*funcs[f])(dim, Nf, NfAux, uOff, uOff_x, u, NULL, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, x, &values[v]); 550 } 551 } else { 552 for (comp = 0; comp < Ncf[f]; ++comp) values[v+comp] = 0.0; 553 } 554 v += Ncf[f]; 555 } 556 } 557 ierr = DMPlexVecRestoreClosure(dm, section, localU, c, NULL, &coefficients);CHKERRQ(ierr); 558 if (dmAux) {ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);CHKERRQ(ierr);} 559 ierr = DMPlexVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr); 560 } 561 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 562 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 563 for (f = 0; f < Nf; f++) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);} 564 ierr = PetscFree2(sp, Ncf);CHKERRQ(ierr); 565 PetscFunctionReturn(0); 566 } 567 568 #undef __FUNCT__ 569 #define __FUNCT__ "DMPlexInsertBoundaryValues_FEM_Internal" 570 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) 571 { 572 PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx); 573 void **ctxs; 574 PetscInt numFields; 575 PetscErrorCode ierr; 576 577 PetscFunctionBegin; 578 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 579 ierr = PetscCalloc2(numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 580 funcs[field] = func; 581 ctxs[field] = ctx; 582 ierr = DMProjectFunctionLabelLocal(dm, time, label, numids, ids, funcs, ctxs, INSERT_BC_VALUES, locX);CHKERRQ(ierr); 583 ierr = PetscFree2(funcs,ctxs);CHKERRQ(ierr); 584 PetscFunctionReturn(0); 585 } 586 587 #undef __FUNCT__ 588 #define __FUNCT__ "DMPlexInsertBoundaryValues_FVM_Internal" 589 /* This ignores numcomps/comps */ 590 static PetscErrorCode DMPlexInsertBoundaryValues_FVM_Internal(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad, 591 PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX) 592 { 593 PetscDS prob; 594 PetscSF sf; 595 DM dmFace, dmCell, dmGrad; 596 const PetscScalar *facegeom, *cellgeom = NULL, *grad; 597 const PetscInt *leaves; 598 PetscScalar *x, *fx; 599 PetscInt dim, nleaves, loc, fStart, fEnd, pdim, i; 600 PetscErrorCode ierr; 601 602 PetscFunctionBegin; 603 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 604 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 605 nleaves = PetscMax(0, nleaves); 606 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 607 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 608 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 609 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 610 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 611 if (cellGeometry) { 612 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 613 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 614 } 615 if (Grad) { 616 PetscFV fv; 617 618 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);CHKERRQ(ierr); 619 ierr = VecGetDM(Grad, &dmGrad);CHKERRQ(ierr); 620 ierr = VecGetArrayRead(Grad, &grad);CHKERRQ(ierr); 621 ierr = PetscFVGetNumComponents(fv, &pdim);CHKERRQ(ierr); 622 ierr = DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 623 } 624 ierr = VecGetArray(locX, &x);CHKERRQ(ierr); 625 for (i = 0; i < numids; ++i) { 626 IS faceIS; 627 const PetscInt *faces; 628 PetscInt numFaces, f; 629 630 ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr); 631 if (!faceIS) continue; /* No points with that id on this process */ 632 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 633 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 634 for (f = 0; f < numFaces; ++f) { 635 const PetscInt face = faces[f], *cells; 636 const PetscFVFaceGeom *fg; 637 638 if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */ 639 ierr = PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);CHKERRQ(ierr); 640 if (loc >= 0) continue; 641 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 642 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 643 if (Grad) { 644 const PetscFVCellGeom *cg; 645 const PetscScalar *cx, *cgrad; 646 PetscScalar *xG; 647 PetscReal dx[3]; 648 PetscInt d; 649 650 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);CHKERRQ(ierr); 651 ierr = DMPlexPointLocalRead(dm, cells[0], x, &cx);CHKERRQ(ierr); 652 ierr = DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);CHKERRQ(ierr); 653 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 654 DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx); 655 for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx); 656 ierr = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);CHKERRQ(ierr); 657 } else { 658 const PetscScalar *xI; 659 PetscScalar *xG; 660 661 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr); 662 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 663 ierr = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);CHKERRQ(ierr); 664 } 665 } 666 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 667 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 668 } 669 ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr); 670 if (Grad) { 671 ierr = DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 672 ierr = VecRestoreArrayRead(Grad, &grad);CHKERRQ(ierr); 673 } 674 if (cellGeometry) { 675 ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 676 } 677 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 678 PetscFunctionReturn(0); 679 } 680 681 #undef __FUNCT__ 682 #define __FUNCT__ "DMPlexInsertBoundaryValues" 683 PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM) 684 { 685 PetscInt numBd, b; 686 PetscErrorCode ierr; 687 688 PetscFunctionBegin; 689 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 690 PetscValidHeaderSpecific(locX, VEC_CLASSID, 2); 691 if (faceGeomFVM) {PetscValidHeaderSpecific(faceGeomFVM, VEC_CLASSID, 4);} 692 if (cellGeomFVM) {PetscValidHeaderSpecific(cellGeomFVM, VEC_CLASSID, 5);} 693 if (gradFVM) {PetscValidHeaderSpecific(gradFVM, VEC_CLASSID, 6);} 694 ierr = DMGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 695 for (b = 0; b < numBd; ++b) { 696 PetscBool isEssential; 697 const char *labelname; 698 DMLabel label; 699 PetscInt field; 700 PetscObject obj; 701 PetscClassId id; 702 void (*func)(); 703 PetscInt numids; 704 const PetscInt *ids; 705 void *ctx; 706 707 ierr = DMGetBoundary(dm, b, &isEssential, NULL, &labelname, &field, NULL, NULL, &func, &numids, &ids, &ctx);CHKERRQ(ierr); 708 if (insertEssential != isEssential) continue; 709 ierr = DMGetLabel(dm, labelname, &label);CHKERRQ(ierr); 710 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 711 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 712 if (id == PETSCFE_CLASSID) { 713 if (!isEssential) continue; /* for FEM, there is no insertion to be done for non-essential boundary conditions */ 714 ierr = DMPlexInsertBoundaryValues_FEM_Internal(dm, time, field, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);CHKERRQ(ierr); 715 } else if (id == PETSCFV_CLASSID) { 716 if (!faceGeomFVM) continue; 717 ierr = DMPlexInsertBoundaryValues_FVM_Internal(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, 718 field, label, numids, ids, (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);CHKERRQ(ierr); 719 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 720 } 721 PetscFunctionReturn(0); 722 } 723 724 #undef __FUNCT__ 725 #define __FUNCT__ "DMComputeL2Diff_Plex" 726 PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 727 { 728 const PetscInt debug = 0; 729 PetscSection section; 730 PetscQuadrature quad; 731 Vec localX; 732 PetscScalar *funcVal, *interpolant; 733 PetscReal *coords, *v0, *J, *invJ, detJ; 734 PetscReal localDiff = 0.0; 735 const PetscReal *quadPoints, *quadWeights; 736 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 737 PetscErrorCode ierr; 738 739 PetscFunctionBegin; 740 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 741 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 742 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 743 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 744 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 745 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 746 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 747 for (field = 0; field < numFields; ++field) { 748 PetscObject obj; 749 PetscClassId id; 750 PetscInt Nc; 751 752 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 753 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 754 if (id == PETSCFE_CLASSID) { 755 PetscFE fe = (PetscFE) obj; 756 757 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 758 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 759 } else if (id == PETSCFV_CLASSID) { 760 PetscFV fv = (PetscFV) obj; 761 762 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 763 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 764 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 765 numComponents += Nc; 766 } 767 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 768 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 769 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 770 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 771 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 772 for (c = cStart; c < cEnd; ++c) { 773 PetscScalar *x = NULL; 774 PetscReal elemDiff = 0.0; 775 776 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 777 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 778 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 779 780 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 781 PetscObject obj; 782 PetscClassId id; 783 void * const ctx = ctxs ? ctxs[field] : NULL; 784 PetscInt Nb, Nc, q, fc; 785 786 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 787 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 788 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 789 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 790 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 791 if (debug) { 792 char title[1024]; 793 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 794 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 795 } 796 for (q = 0; q < numQuadPoints; ++q) { 797 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 798 ierr = (*funcs[field])(dim, time, coords, Nc, funcVal, ctx);CHKERRQ(ierr); 799 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 800 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 801 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 802 for (fc = 0; fc < Nc; ++fc) { 803 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);} 804 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 805 } 806 } 807 fieldOffset += Nb*Nc; 808 } 809 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 810 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 811 localDiff += elemDiff; 812 } 813 ierr = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 814 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 815 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 816 *diff = PetscSqrtReal(*diff); 817 PetscFunctionReturn(0); 818 } 819 820 #undef __FUNCT__ 821 #define __FUNCT__ "DMComputeL2GradientDiff_Plex" 822 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) 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, time, 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__ "DMComputeL2FieldDiff_Plex" 924 PetscErrorCode DMComputeL2FieldDiff_Plex(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 = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 943 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 944 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 945 for (field = 0; field < numFields; ++field) { 946 PetscObject obj; 947 PetscClassId id; 948 PetscInt Nc; 949 950 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 951 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 952 if (id == PETSCFE_CLASSID) { 953 PetscFE fe = (PetscFE) obj; 954 955 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 956 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 957 } else if (id == PETSCFV_CLASSID) { 958 PetscFV fv = (PetscFV) obj; 959 960 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 961 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 962 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 963 numComponents += Nc; 964 } 965 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);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 = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 1055 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1056 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1057 for (field = 0; field < numFields; ++field) { 1058 PetscObject obj; 1059 PetscClassId id; 1060 PetscInt Nc; 1061 1062 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1063 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1064 if (id == PETSCFE_CLASSID) { 1065 PetscFE fe = (PetscFE) obj; 1066 1067 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1068 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1069 } else if (id == PETSCFV_CLASSID) { 1070 PetscFV fv = (PetscFV) obj; 1071 1072 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1073 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1074 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1075 numComponents += Nc; 1076 } 1077 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);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 = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 1157 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1158 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1159 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1160 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1161 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1162 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1163 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1164 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1165 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1166 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1167 numCells = cEnd - cStart; 1168 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1169 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1170 if (dmAux) { 1171 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1172 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1173 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1174 } 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 Mat preallocator; 1406 PetscScalar *vals; 1407 PetscInt *cellCIndices, *cellFIndices; 1408 PetscInt locRows, locCols, cell; 1409 1410 ierr = MatGetLocalSize(In, &locRows, &locCols);CHKERRQ(ierr); 1411 ierr = MatCreate(PetscObjectComm((PetscObject) In), &preallocator);CHKERRQ(ierr); 1412 ierr = MatSetType(preallocator, MATPREALLOCATOR);CHKERRQ(ierr); 1413 ierr = MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1414 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 1415 ierr = PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1416 for (cell = cStart; cell < cEnd; ++cell) { 1417 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1418 ierr = MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);CHKERRQ(ierr); 1419 } 1420 ierr = PetscFree3(vals,cellCIndices,cellFIndices);CHKERRQ(ierr); 1421 ierr = MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1422 ierr = MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1423 ierr = MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);CHKERRQ(ierr); 1424 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 1425 } 1426 /* Fill matrix */ 1427 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1428 for (c = cStart; c < cEnd; ++c) { 1429 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1430 } 1431 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1432 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1433 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1434 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1435 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1436 if (mesh->printFEM) { 1437 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1438 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1439 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1440 } 1441 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1442 PetscFunctionReturn(0); 1443 } 1444 1445 #undef __FUNCT__ 1446 #define __FUNCT__ "DMPlexComputeInterpolatorGeneral" 1447 /*@ 1448 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 1449 1450 Input Parameters: 1451 + dmf - The fine mesh 1452 . dmc - The coarse mesh 1453 - user - The user context 1454 1455 Output Parameter: 1456 . In - The interpolation matrix 1457 1458 Level: developer 1459 1460 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1461 @*/ 1462 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 1463 { 1464 DM_Plex *mesh = (DM_Plex *) dmf->data; 1465 const char *name = "Interpolator"; 1466 PetscDS prob; 1467 PetscSection fsection, csection, globalFSection, globalCSection; 1468 PetscHashJK ht; 1469 PetscLayout rLayout; 1470 PetscInt *dnz, *onz; 1471 PetscInt locRows, rStart, rEnd; 1472 PetscReal *x, *v0, *J, *invJ, detJ; 1473 PetscReal *v0c, *Jc, *invJc, detJc; 1474 PetscScalar *elemMat; 1475 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1476 PetscErrorCode ierr; 1477 1478 PetscFunctionBegin; 1479 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1480 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1481 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1482 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1483 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1484 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1485 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1486 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1487 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1488 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1489 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1490 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1491 ierr = PetscMalloc1(totDim*totDim, &elemMat);CHKERRQ(ierr); 1492 1493 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1494 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1495 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1496 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1497 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1498 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1499 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1500 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1501 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1502 for (field = 0; field < Nf; ++field) { 1503 PetscObject obj; 1504 PetscClassId id; 1505 PetscDualSpace Q = NULL; 1506 PetscQuadrature f; 1507 const PetscReal *qpoints; 1508 PetscInt Nc, Np, fpdim, i, d; 1509 1510 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1511 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1512 if (id == PETSCFE_CLASSID) { 1513 PetscFE fe = (PetscFE) obj; 1514 1515 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1516 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1517 } else if (id == PETSCFV_CLASSID) { 1518 PetscFV fv = (PetscFV) obj; 1519 1520 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1521 Nc = 1; 1522 } 1523 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1524 /* For each fine grid cell */ 1525 for (cell = cStart; cell < cEnd; ++cell) { 1526 PetscInt *findices, *cindices; 1527 PetscInt numFIndices, numCIndices; 1528 1529 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1530 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1531 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1532 for (i = 0; i < fpdim; ++i) { 1533 Vec pointVec; 1534 PetscScalar *pV; 1535 IS coarseCellIS; 1536 const PetscInt *coarseCells; 1537 PetscInt numCoarseCells, q, r, c; 1538 1539 /* Get points from the dual basis functional quadrature */ 1540 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1541 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1542 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1543 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1544 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1545 for (q = 0; q < Np; ++q) { 1546 /* Transform point to real space */ 1547 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1548 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1549 } 1550 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1551 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1552 ierr = DMLocatePoints(dmc, pointVec, &coarseCellIS);CHKERRQ(ierr); 1553 ierr = ISViewFromOptions(coarseCellIS, NULL, "-interp_is_view");CHKERRQ(ierr); 1554 /* Update preallocation info */ 1555 ierr = ISGetLocalSize(coarseCellIS, &numCoarseCells);CHKERRQ(ierr); 1556 ierr = ISGetIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1557 for (r = 0; r < Nc; ++r) { 1558 PetscHashJKKey key; 1559 PetscHashJKIter missing, iter; 1560 1561 key.j = findices[i*Nc+r]; 1562 if (key.j < 0) continue; 1563 /* Get indices for coarse elements */ 1564 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1565 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1566 for (c = 0; c < numCIndices; ++c) { 1567 key.k = cindices[c]; 1568 if (key.k < 0) continue; 1569 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1570 if (missing) { 1571 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1572 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1573 else ++onz[key.j-rStart]; 1574 } 1575 } 1576 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1577 } 1578 } 1579 ierr = ISRestoreIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1580 ierr = ISDestroy(&coarseCellIS);CHKERRQ(ierr); 1581 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1582 } 1583 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1584 } 1585 } 1586 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1587 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1588 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1589 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1590 for (field = 0; field < Nf; ++field) { 1591 PetscObject obj; 1592 PetscClassId id; 1593 PetscDualSpace Q = NULL; 1594 PetscQuadrature f; 1595 const PetscReal *qpoints, *qweights; 1596 PetscInt Nc, Np, fpdim, i, d; 1597 1598 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1599 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1600 if (id == PETSCFE_CLASSID) { 1601 PetscFE fe = (PetscFE) obj; 1602 1603 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1604 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1605 } else if (id == PETSCFV_CLASSID) { 1606 PetscFV fv = (PetscFV) obj; 1607 1608 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1609 Nc = 1; 1610 } 1611 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1612 /* For each fine grid cell */ 1613 for (cell = cStart; cell < cEnd; ++cell) { 1614 PetscInt *findices, *cindices; 1615 PetscInt numFIndices, numCIndices; 1616 1617 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1618 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1619 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1620 for (i = 0; i < fpdim; ++i) { 1621 Vec pointVec; 1622 PetscScalar *pV; 1623 IS coarseCellIS; 1624 const PetscInt *coarseCells; 1625 PetscInt numCoarseCells, cpdim, q, c, j; 1626 1627 /* Get points from the dual basis functional quadrature */ 1628 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1629 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, &qweights);CHKERRQ(ierr); 1630 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1631 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1632 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1633 for (q = 0; q < Np; ++q) { 1634 /* Transform point to real space */ 1635 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1636 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1637 } 1638 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1639 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1640 ierr = DMLocatePoints(dmc, pointVec, &coarseCellIS);CHKERRQ(ierr); 1641 /* Update preallocation info */ 1642 ierr = ISGetLocalSize(coarseCellIS, &numCoarseCells);CHKERRQ(ierr); 1643 ierr = ISGetIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1644 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1645 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1646 PetscReal pVReal[3]; 1647 1648 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1649 /* Transform points from real space to coarse reference space */ 1650 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell], NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 1651 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 1652 CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x); 1653 1654 if (id == PETSCFE_CLASSID) { 1655 PetscFE fe = (PetscFE) obj; 1656 PetscReal *B; 1657 1658 /* Evaluate coarse basis on contained point */ 1659 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1660 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 1661 /* Get elemMat entries by multiplying by weight */ 1662 for (j = 0; j < cpdim; ++j) { 1663 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = B[j*Nc + c]*qweights[ccell]; 1664 } 1665 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1666 } else { 1667 cpdim = 1; 1668 for (j = 0; j < cpdim; ++j) { 1669 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = 1.0*qweights[ccell]; 1670 } 1671 } 1672 /* Update interpolator */ 1673 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, Nc, numCIndices, elemMat);CHKERRQ(ierr);} 1674 ierr = MatSetValues(In, Nc, &findices[i*Nc], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1675 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1676 } 1677 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1678 ierr = ISRestoreIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1679 ierr = ISDestroy(&coarseCellIS);CHKERRQ(ierr); 1680 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1681 } 1682 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1683 } 1684 } 1685 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 1686 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 1687 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1688 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1689 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1690 PetscFunctionReturn(0); 1691 } 1692 1693 #undef __FUNCT__ 1694 #define __FUNCT__ "DMPlexComputeInjectorFEM" 1695 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 1696 { 1697 PetscDS prob; 1698 PetscFE *feRef; 1699 PetscFV *fvRef; 1700 Vec fv, cv; 1701 IS fis, cis; 1702 PetscSection fsection, fglobalSection, csection, cglobalSection; 1703 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 1704 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m; 1705 PetscErrorCode ierr; 1706 1707 PetscFunctionBegin; 1708 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1709 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1710 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1711 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1712 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1713 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1714 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1715 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1716 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1717 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1718 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1719 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1720 for (f = 0; f < Nf; ++f) { 1721 PetscObject obj; 1722 PetscClassId id; 1723 PetscInt fNb = 0, Nc = 0; 1724 1725 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1726 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1727 if (id == PETSCFE_CLASSID) { 1728 PetscFE fe = (PetscFE) obj; 1729 1730 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1731 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 1732 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1733 } else if (id == PETSCFV_CLASSID) { 1734 PetscFV fv = (PetscFV) obj; 1735 PetscDualSpace Q; 1736 1737 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1738 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1739 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 1740 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1741 } 1742 fTotDim += fNb*Nc; 1743 } 1744 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1745 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 1746 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 1747 PetscFE feC; 1748 PetscFV fvC; 1749 PetscDualSpace QF, QC; 1750 PetscInt NcF, NcC, fpdim, cpdim; 1751 1752 if (feRef[field]) { 1753 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 1754 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 1755 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 1756 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 1757 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1758 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 1759 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1760 } else { 1761 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 1762 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 1763 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 1764 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 1765 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1766 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 1767 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1768 } 1769 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); 1770 for (c = 0; c < cpdim; ++c) { 1771 PetscQuadrature cfunc; 1772 const PetscReal *cqpoints; 1773 PetscInt NpC; 1774 PetscBool found = PETSC_FALSE; 1775 1776 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 1777 ierr = PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 1778 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 1779 for (f = 0; f < fpdim; ++f) { 1780 PetscQuadrature ffunc; 1781 const PetscReal *fqpoints; 1782 PetscReal sum = 0.0; 1783 PetscInt NpF, comp; 1784 1785 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 1786 ierr = PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 1787 if (NpC != NpF) continue; 1788 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 1789 if (sum > 1.0e-9) continue; 1790 for (comp = 0; comp < NcC; ++comp) { 1791 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 1792 } 1793 found = PETSC_TRUE; 1794 break; 1795 } 1796 if (!found) { 1797 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 1798 if (fvRef[field]) { 1799 PetscInt comp; 1800 for (comp = 0; comp < NcC; ++comp) { 1801 cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp; 1802 } 1803 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 1804 } 1805 } 1806 offsetC += cpdim*NcC; 1807 offsetF += fpdim*NcF; 1808 } 1809 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 1810 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1811 1812 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 1813 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 1814 ierr = VecGetOwnershipRange(cv, &startC, &endC);CHKERRQ(ierr); 1815 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 1816 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 1817 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 1818 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 1819 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 1820 for (c = cStart; c < cEnd; ++c) { 1821 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 1822 for (d = 0; d < cTotDim; ++d) { 1823 if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue; 1824 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]]); 1825 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 1826 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 1827 } 1828 } 1829 ierr = PetscFree(cmap);CHKERRQ(ierr); 1830 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 1831 1832 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 1833 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 1834 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 1835 ierr = ISDestroy(&cis);CHKERRQ(ierr); 1836 ierr = ISDestroy(&fis);CHKERRQ(ierr); 1837 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 1838 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 1839 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1840 PetscFunctionReturn(0); 1841 } 1842