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