1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petscsf.h> 3 4 #include <petsc/private/petscfeimpl.h> 5 #include <petsc/private/petscfvimpl.h> 6 7 #undef __FUNCT__ 8 #define __FUNCT__ "DMPlexGetScale" 9 PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale) 10 { 11 DM_Plex *mesh = (DM_Plex*) dm->data; 12 13 PetscFunctionBegin; 14 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 15 PetscValidPointer(scale, 3); 16 *scale = mesh->scale[unit]; 17 PetscFunctionReturn(0); 18 } 19 20 #undef __FUNCT__ 21 #define __FUNCT__ "DMPlexSetScale" 22 PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale) 23 { 24 DM_Plex *mesh = (DM_Plex*) dm->data; 25 26 PetscFunctionBegin; 27 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 28 mesh->scale[unit] = scale; 29 PetscFunctionReturn(0); 30 } 31 32 PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k) 33 { 34 switch (i) { 35 case 0: 36 switch (j) { 37 case 0: return 0; 38 case 1: 39 switch (k) { 40 case 0: return 0; 41 case 1: return 0; 42 case 2: return 1; 43 } 44 case 2: 45 switch (k) { 46 case 0: return 0; 47 case 1: return -1; 48 case 2: return 0; 49 } 50 } 51 case 1: 52 switch (j) { 53 case 0: 54 switch (k) { 55 case 0: return 0; 56 case 1: return 0; 57 case 2: return -1; 58 } 59 case 1: return 0; 60 case 2: 61 switch (k) { 62 case 0: return 1; 63 case 1: return 0; 64 case 2: return 0; 65 } 66 } 67 case 2: 68 switch (j) { 69 case 0: 70 switch (k) { 71 case 0: return 0; 72 case 1: return 1; 73 case 2: return 0; 74 } 75 case 1: 76 switch (k) { 77 case 0: return -1; 78 case 1: return 0; 79 case 2: return 0; 80 } 81 case 2: return 0; 82 } 83 } 84 return 0; 85 } 86 87 #undef __FUNCT__ 88 #define __FUNCT__ "DMPlexProjectRigidBody" 89 static PetscErrorCode DMPlexProjectRigidBody(PetscInt dim, PetscReal t, const PetscReal X[], PetscInt Nf, PetscScalar *mode, void *ctx) 90 { 91 PetscInt *ctxInt = (PetscInt *) ctx; 92 PetscInt dim2 = ctxInt[0]; 93 PetscInt d = ctxInt[1]; 94 PetscInt i, j, k = dim > 2 ? d - dim : d; 95 96 PetscFunctionBegin; 97 if (dim != dim2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input dimension %d does not match context dimension %d", dim, dim2); 98 for (i = 0; i < dim; i++) mode[i] = 0.; 99 if (d < dim) { 100 mode[d] = 1.; 101 } else { 102 for (i = 0; i < dim; i++) { 103 for (j = 0; j < dim; j++) { 104 mode[j] += epsilon(i, j, k)*X[i]; 105 } 106 } 107 } 108 PetscFunctionReturn(0); 109 } 110 111 #undef __FUNCT__ 112 #define __FUNCT__ "DMPlexCreateRigidBody" 113 /*@C 114 DMPlexCreateRigidBody - for the default global section, create rigid body modes from coordinates 115 116 Collective on DM 117 118 Input Arguments: 119 . dm - the DM 120 121 Output Argument: 122 . sp - the null space 123 124 Note: This is necessary to take account of Dirichlet conditions on the displacements 125 126 Level: advanced 127 128 .seealso: MatNullSpaceCreate() 129 @*/ 130 PetscErrorCode DMPlexCreateRigidBody(DM dm, MatNullSpace *sp) 131 { 132 MPI_Comm comm; 133 Vec mode[6]; 134 PetscSection section, globalSection; 135 PetscInt dim, dimEmbed, n, m, d, i, j; 136 PetscErrorCode ierr; 137 138 PetscFunctionBegin; 139 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 140 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 141 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 142 if (dim == 1) { 143 ierr = MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);CHKERRQ(ierr); 144 PetscFunctionReturn(0); 145 } 146 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 147 ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr); 148 ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr); 149 m = (dim*(dim+1))/2; 150 ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr); 151 ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr); 152 ierr = VecSetUp(mode[0]);CHKERRQ(ierr); 153 for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);} 154 for (d = 0; d < m; d++) { 155 PetscInt ctx[2]; 156 void *voidctx = (void *) (&ctx[0]); 157 PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody; 158 159 ctx[0] = dimEmbed; 160 ctx[1] = d; 161 ierr = DMProjectFunction(dm, 0.0, &func, &voidctx, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 162 } 163 for (i = 0; i < dim; ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);} 164 /* Orthonormalize system */ 165 for (i = dim; i < m; ++i) { 166 PetscScalar dots[6]; 167 168 ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr); 169 for (j = 0; j < i; ++j) dots[j] *= -1.0; 170 ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr); 171 ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr); 172 } 173 ierr = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);CHKERRQ(ierr); 174 for (i = 0; i< m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);} 175 PetscFunctionReturn(0); 176 } 177 178 #undef __FUNCT__ 179 #define __FUNCT__ "DMPlexSetMaxProjectionHeight" 180 /*@ 181 DMPlexSetMaxProjectionHeight - In DMPlexProjectXXXLocal() functions, the projected values of a basis function's dofs 182 are computed by associating the basis function with one of the mesh points in its transitively-closed support, and 183 evaluating the dual space basis of that point. A basis function is associated with the point in its 184 transitively-closed support whose mesh height is highest (w.r.t. DAG height), but not greater than the maximum 185 projection height, which is set with this function. By default, the maximum projection height is zero, which means 186 that only mesh cells are used to project basis functions. A height of one, for example, evaluates a cell-interior 187 basis functions using its cells dual space basis, but all other basis functions with the dual space basis of a face. 188 189 Input Parameters: 190 + dm - the DMPlex object 191 - height - the maximum projection height >= 0 192 193 Level: advanced 194 195 .seealso: DMPlexGetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal() 196 @*/ 197 PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height) 198 { 199 DM_Plex *plex = (DM_Plex *) dm->data; 200 201 PetscFunctionBegin; 202 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 203 plex->maxProjectionHeight = height; 204 PetscFunctionReturn(0); 205 } 206 207 #undef __FUNCT__ 208 #define __FUNCT__ "DMPlexGetMaxProjectionHeight" 209 /*@ 210 DMPlexGetMaxProjectionHeight - Get the maximum height (w.r.t. DAG) of mesh points used to evaluate dual bases in 211 DMPlexProjectXXXLocal() functions. 212 213 Input Parameters: 214 . dm - the DMPlex object 215 216 Output Parameters: 217 . height - the maximum projection height 218 219 Level: intermediate 220 221 .seealso: DMPlexSetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal() 222 @*/ 223 PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height) 224 { 225 DM_Plex *plex = (DM_Plex *) dm->data; 226 227 PetscFunctionBegin; 228 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 229 *height = plex->maxProjectionHeight; 230 PetscFunctionReturn(0); 231 } 232 233 #undef __FUNCT__ 234 #define __FUNCT__ "DMProjectFunctionLabelLocal_Plex" 235 PetscErrorCode DMProjectFunctionLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 236 { 237 PetscDualSpace *sp, *cellsp; 238 PetscInt *numComp; 239 PetscSection section; 240 PetscScalar *values; 241 PetscBool *fieldActive; 242 PetscInt numFields, dim, dimEmbed, spDim, totDim = 0, numValues, pStart, pEnd, cStart, cEnd, cEndInterior, f, d, v, i, comp, maxHeight, h; 243 PetscErrorCode ierr; 244 245 PetscFunctionBegin; 246 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 247 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 248 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 249 if (cEnd <= cStart) PetscFunctionReturn(0); 250 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 251 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 252 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 253 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 254 ierr = PetscMalloc2(numFields,&sp,numFields,&numComp);CHKERRQ(ierr); 255 ierr = DMPlexGetMaxProjectionHeight(dm,&maxHeight);CHKERRQ(ierr); 256 if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"maximum projection height %d not in [0, %d)\n", maxHeight,dim);} 257 if (maxHeight > 0) {ierr = PetscMalloc1(numFields,&cellsp);CHKERRQ(ierr);} 258 else {cellsp = sp;} 259 for (h = 0; h <= maxHeight; h++) { 260 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 261 if (!h) {pStart = cStart; pEnd = cEnd;} 262 if (pEnd <= pStart) continue; 263 totDim = 0; 264 for (f = 0; f < numFields; ++f) { 265 PetscObject obj; 266 PetscClassId id; 267 268 ierr = DMGetField(dm, f, &obj);CHKERRQ(ierr); 269 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 270 if (id == PETSCFE_CLASSID) { 271 PetscFE fe = (PetscFE) obj; 272 273 ierr = PetscFEGetNumComponents(fe, &numComp[f]);CHKERRQ(ierr); 274 if (!h) { 275 ierr = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr); 276 sp[f] = cellsp[f]; 277 ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr); 278 } else { 279 ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);CHKERRQ(ierr); 280 if (!sp[f]) continue; 281 } 282 } else if (id == PETSCFV_CLASSID) { 283 PetscFV fv = (PetscFV) obj; 284 285 ierr = PetscFVGetNumComponents(fv, &numComp[f]);CHKERRQ(ierr); 286 ierr = PetscFVGetDualSpace(fv, &sp[f]);CHKERRQ(ierr); 287 ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr); 288 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 289 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 290 totDim += spDim*numComp[f]; 291 } 292 ierr = DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);CHKERRQ(ierr); 293 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim); 294 if (!totDim) continue; 295 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 296 ierr = DMGetWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);CHKERRQ(ierr); 297 for (f = 0; f < numFields; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE; 298 for (i = 0; i < numIds; ++i) { 299 IS pointIS; 300 const PetscInt *points; 301 PetscInt n, p; 302 303 ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr); 304 ierr = ISGetLocalSize(pointIS, &n);CHKERRQ(ierr); 305 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 306 for (p = 0; p < n; ++p) { 307 const PetscInt point = points[p]; 308 PetscFECellGeom geom; 309 310 if ((point < pStart) || (point >= pEnd)) continue; 311 ierr = DMPlexComputeCellGeometryFEM(dm, point, NULL, geom.v0, geom.J, NULL, &geom.detJ);CHKERRQ(ierr); 312 geom.dim = dim - h; 313 geom.dimEmbed = dimEmbed; 314 for (f = 0, v = 0; f < numFields; ++f) { 315 void * const ctx = ctxs ? ctxs[f] : NULL; 316 317 if (!sp[f]) continue; 318 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 319 for (d = 0; d < spDim; ++d) { 320 if (funcs[f]) { 321 ierr = PetscDualSpaceApply(sp[f], d, time, &geom, numComp[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr); 322 } else { 323 for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0; 324 } 325 v += numComp[f]; 326 } 327 } 328 ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, values, mode);CHKERRQ(ierr); 329 } 330 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 331 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 332 } 333 if (h) { 334 for (f = 0; f < numFields; ++f) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);} 335 } 336 } 337 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 338 ierr = DMRestoreWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);CHKERRQ(ierr); 339 for (f = 0; f < numFields; ++f) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);} 340 ierr = PetscFree2(sp, numComp);CHKERRQ(ierr); 341 if (maxHeight > 0) { 342 ierr = PetscFree(cellsp);CHKERRQ(ierr); 343 } 344 PetscFunctionReturn(0); 345 } 346 347 #undef __FUNCT__ 348 #define __FUNCT__ "DMProjectFunctionLocal_Plex" 349 PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 350 { 351 PetscDualSpace *sp, *cellsp; 352 PetscInt *numComp; 353 PetscSection section; 354 PetscScalar *values; 355 PetscInt Nf, dim, dimEmbed, spDim, totDim = 0, numValues, pStart, pEnd, p, cStart, cEnd, cEndInterior, f, d, v, comp, h, maxHeight; 356 PetscErrorCode ierr; 357 358 PetscFunctionBegin; 359 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 360 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 361 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 362 if (cEnd <= cStart) PetscFunctionReturn(0); 363 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 364 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 365 ierr = PetscMalloc2(Nf, &sp, Nf, &numComp);CHKERRQ(ierr); 366 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 367 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 368 ierr = DMPlexGetMaxProjectionHeight(dm,&maxHeight);CHKERRQ(ierr); 369 if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"maximum projection height %d not in [0, %d)\n", maxHeight,dim);} 370 if (maxHeight > 0) { 371 ierr = PetscMalloc1(Nf,&cellsp);CHKERRQ(ierr); 372 } 373 else { 374 cellsp = sp; 375 } 376 for (h = 0; h <= maxHeight; h++) { 377 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 378 if (!h) {pStart = cStart; pEnd = cEnd;} 379 if (pEnd <= pStart) continue; 380 totDim = 0; 381 for (f = 0; f < Nf; ++f) { 382 PetscObject obj; 383 PetscClassId id; 384 385 ierr = DMGetField(dm, f, &obj);CHKERRQ(ierr); 386 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 387 if (id == PETSCFE_CLASSID) { 388 PetscFE fe = (PetscFE) obj; 389 390 ierr = PetscFEGetNumComponents(fe, &numComp[f]);CHKERRQ(ierr); 391 if (!h) { 392 ierr = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr); 393 sp[f] = cellsp[f]; 394 ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr); 395 } 396 else { 397 ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);CHKERRQ(ierr); 398 if (!sp[f]) { 399 continue; 400 } 401 } 402 } else if (id == PETSCFV_CLASSID) { 403 PetscFV fv = (PetscFV) obj; 404 405 if (!h) { 406 ierr = PetscFVGetNumComponents(fv, &numComp[f]);CHKERRQ(ierr); 407 ierr = PetscFVGetDualSpace(fv, &cellsp[f]);CHKERRQ(ierr); 408 ierr = PetscObjectReference((PetscObject) cellsp[f]);CHKERRQ(ierr); 409 sp[f] = cellsp[f]; 410 } 411 else { 412 sp[f] = NULL; 413 continue; /* FV doesn't have fields that can be evaluated at faces, corners, etc. */ 414 } 415 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 416 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 417 totDim += spDim*numComp[f]; 418 } 419 ierr = DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);CHKERRQ(ierr); 420 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim); 421 if (!totDim) continue; 422 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 423 for (p = pStart; p < pEnd; ++p) { 424 PetscFECellGeom geom; 425 426 ierr = DMPlexComputeCellGeometryFEM(dm, p, NULL, geom.v0, geom.J, NULL, &geom.detJ);CHKERRQ(ierr); 427 geom.dim = dim - h; 428 geom.dimEmbed = dimEmbed; 429 for (f = 0, v = 0; f < Nf; ++f) { 430 void * const ctx = ctxs ? ctxs[f] : NULL; 431 432 if (!sp[f]) continue; 433 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 434 for (d = 0; d < spDim; ++d) { 435 if (funcs[f]) { 436 ierr = PetscDualSpaceApply(sp[f], d, time, &geom, numComp[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr); 437 } else { 438 for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0; 439 } 440 v += numComp[f]; 441 } 442 } 443 ierr = DMPlexVecSetClosure(dm, section, localX, p, values, mode);CHKERRQ(ierr); 444 } 445 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 446 if (h || !maxHeight) { 447 for (f = 0; f < Nf; f++) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);} 448 } 449 } 450 ierr = PetscFree2(sp, numComp);CHKERRQ(ierr); 451 if (maxHeight > 0) { 452 for (f = 0; f < Nf; f++) {ierr = PetscDualSpaceDestroy(&cellsp[f]);CHKERRQ(ierr);} 453 ierr = PetscFree(cellsp);CHKERRQ(ierr); 454 } 455 PetscFunctionReturn(0); 456 } 457 458 #undef __FUNCT__ 459 #define __FUNCT__ "DMPlexProjectFieldLocal" 460 PetscErrorCode DMPlexProjectFieldLocal(DM dm, Vec localU, void (**funcs)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal [], PetscScalar []), InsertMode mode, Vec localX) 461 { 462 DM dmAux; 463 PetscDS prob, probAux = NULL; 464 Vec A; 465 PetscSection section, sectionAux = NULL; 466 PetscDualSpace *sp; 467 PetscInt *Ncf; 468 PetscScalar *values, *u, *u_x, *a, *a_x; 469 PetscReal *x, *v0, *J, *invJ, detJ; 470 PetscInt Nf, dim, spDim, totDim, numValues, cStart, cEnd, cEndInterior, c, f, d, v, comp, maxHeight; 471 PetscErrorCode ierr; 472 473 PetscFunctionBegin; 474 ierr = DMPlexGetMaxProjectionHeight(dm,&maxHeight);CHKERRQ(ierr); 475 if (maxHeight > 0) {SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Field projection for height > 0 not supported yet");} 476 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 477 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 478 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 479 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 480 ierr = PetscMalloc2(Nf, &sp, Nf, &Ncf);CHKERRQ(ierr); 481 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 482 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 483 ierr = PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);CHKERRQ(ierr); 484 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 485 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 486 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 487 if (dmAux) { 488 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 489 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 490 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);CHKERRQ(ierr); 491 } 492 ierr = DMPlexInsertBoundaryValues(dm, localU, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 493 ierr = DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr); 494 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim); 495 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 496 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 497 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 498 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 499 for (c = cStart; c < cEnd; ++c) { 500 PetscScalar *coefficients = NULL, *coefficientsAux = NULL; 501 502 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 503 ierr = DMPlexVecGetClosure(dm, section, localU, c, NULL, &coefficients);CHKERRQ(ierr); 504 if (dmAux) {ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);CHKERRQ(ierr);} 505 for (f = 0, v = 0; f < Nf; ++f) { 506 PetscObject obj; 507 PetscClassId id; 508 509 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 510 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 511 if (id == PETSCFE_CLASSID) { 512 PetscFE fe = (PetscFE) obj; 513 514 ierr = PetscFEGetDualSpace(fe, &sp[f]);CHKERRQ(ierr); 515 ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr); 516 ierr = PetscFEGetNumComponents(fe, &Ncf[f]);CHKERRQ(ierr); 517 } else if (id == PETSCFV_CLASSID) { 518 PetscFV fv = (PetscFV) obj; 519 520 ierr = PetscFVGetNumComponents(fv, &Ncf[f]);CHKERRQ(ierr); 521 ierr = PetscFVGetDualSpace(fv, &sp[f]);CHKERRQ(ierr); 522 ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr); 523 } 524 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 525 for (d = 0; d < spDim; ++d) { 526 PetscQuadrature quad; 527 const PetscReal *points, *weights; 528 PetscInt numPoints, q; 529 530 if (funcs[f]) { 531 ierr = PetscDualSpaceGetFunctional(sp[f], d, &quad);CHKERRQ(ierr); 532 ierr = PetscQuadratureGetData(quad, NULL, &numPoints, &points, &weights);CHKERRQ(ierr); 533 for (q = 0; q < numPoints; ++q) { 534 CoordinatesRefToReal(dim, dim, v0, J, &points[q*dim], x); 535 ierr = EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, coefficients, NULL, u, u_x, NULL);CHKERRQ(ierr); 536 ierr = EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, coefficientsAux, NULL, a, a_x, NULL);CHKERRQ(ierr); 537 (*funcs[f])(u, NULL, u_x, a, NULL, a_x, x, &values[v]); 538 } 539 } else { 540 for (comp = 0; comp < Ncf[f]; ++comp) values[v+comp] = 0.0; 541 } 542 v += Ncf[f]; 543 } 544 } 545 ierr = DMPlexVecRestoreClosure(dm, section, localU, c, NULL, &coefficients);CHKERRQ(ierr); 546 if (dmAux) {ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);CHKERRQ(ierr);} 547 ierr = DMPlexVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr); 548 } 549 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 550 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 551 for (f = 0; f < Nf; f++) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);} 552 ierr = PetscFree2(sp, Ncf);CHKERRQ(ierr); 553 PetscFunctionReturn(0); 554 } 555 556 #undef __FUNCT__ 557 #define __FUNCT__ "DMPlexInsertBoundaryValues_FEM_Internal" 558 static PetscErrorCode DMPlexInsertBoundaryValues_FEM_Internal(DM dm, PetscReal time, PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void *ctx, Vec locX) 559 { 560 PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx); 561 void **ctxs; 562 PetscInt numFields; 563 PetscErrorCode ierr; 564 565 PetscFunctionBegin; 566 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 567 ierr = PetscCalloc2(numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 568 funcs[field] = func; 569 ctxs[field] = ctx; 570 ierr = DMProjectFunctionLabelLocal(dm, time, label, numids, ids, funcs, ctxs, INSERT_BC_VALUES, locX);CHKERRQ(ierr); 571 ierr = PetscFree2(funcs,ctxs);CHKERRQ(ierr); 572 PetscFunctionReturn(0); 573 } 574 575 #undef __FUNCT__ 576 #define __FUNCT__ "DMPlexInsertBoundaryValues_FVM_Internal" 577 /* This ignores numcomps/comps */ 578 static PetscErrorCode DMPlexInsertBoundaryValues_FVM_Internal(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad, 579 PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX) 580 { 581 PetscDS prob; 582 PetscSF sf; 583 DM dmFace, dmCell, dmGrad; 584 const PetscScalar *facegeom, *cellgeom, *grad; 585 const PetscInt *leaves; 586 PetscScalar *x, *fx; 587 PetscInt dim, nleaves, loc, fStart, fEnd, pdim, i; 588 PetscErrorCode ierr; 589 590 PetscFunctionBegin; 591 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 592 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 593 nleaves = PetscMax(0, nleaves); 594 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 595 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 596 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 597 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 598 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 599 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 600 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 601 if (Grad) { 602 PetscFV fv; 603 604 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);CHKERRQ(ierr); 605 ierr = VecGetDM(Grad, &dmGrad);CHKERRQ(ierr); 606 ierr = VecGetArrayRead(Grad, &grad);CHKERRQ(ierr); 607 ierr = PetscFVGetNumComponents(fv, &pdim);CHKERRQ(ierr); 608 ierr = DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 609 } 610 ierr = VecGetArray(locX, &x);CHKERRQ(ierr); 611 for (i = 0; i < numids; ++i) { 612 IS faceIS; 613 const PetscInt *faces; 614 PetscInt numFaces, f; 615 616 ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr); 617 if (!faceIS) continue; /* No points with that id on this process */ 618 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 619 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 620 for (f = 0; f < numFaces; ++f) { 621 const PetscInt face = faces[f], *cells; 622 const PetscFVFaceGeom *fg; 623 624 if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */ 625 ierr = PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);CHKERRQ(ierr); 626 if (loc >= 0) continue; 627 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 628 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 629 if (Grad) { 630 const PetscFVCellGeom *cg; 631 const PetscScalar *cx, *cgrad; 632 PetscScalar *xG; 633 PetscReal dx[3]; 634 PetscInt d; 635 636 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);CHKERRQ(ierr); 637 ierr = DMPlexPointLocalRead(dm, cells[0], x, &cx);CHKERRQ(ierr); 638 ierr = DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);CHKERRQ(ierr); 639 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 640 DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx); 641 for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx); 642 ierr = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);CHKERRQ(ierr); 643 } else { 644 const PetscScalar *xI; 645 PetscScalar *xG; 646 647 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr); 648 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 649 ierr = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);CHKERRQ(ierr); 650 } 651 } 652 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 653 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 654 } 655 ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr); 656 if (Grad) { 657 ierr = DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 658 ierr = VecRestoreArrayRead(Grad, &grad);CHKERRQ(ierr); 659 } 660 ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 661 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 662 PetscFunctionReturn(0); 663 } 664 665 #undef __FUNCT__ 666 #define __FUNCT__ "DMPlexInsertBoundaryValues" 667 PetscErrorCode DMPlexInsertBoundaryValues(DM dm, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM) 668 { 669 PetscInt numBd, b; 670 PetscErrorCode ierr; 671 672 PetscFunctionBegin; 673 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 674 PetscValidHeaderSpecific(locX, VEC_CLASSID, 2); 675 if (faceGeomFVM) {PetscValidHeaderSpecific(faceGeomFVM, VEC_CLASSID, 4);} 676 if (cellGeomFVM) {PetscValidHeaderSpecific(cellGeomFVM, VEC_CLASSID, 5);} 677 if (gradFVM) {PetscValidHeaderSpecific(gradFVM, VEC_CLASSID, 6);} 678 ierr = DMGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 679 for (b = 0; b < numBd; ++b) { 680 PetscBool isEssential; 681 const char *labelname; 682 DMLabel label; 683 PetscInt field; 684 PetscObject obj; 685 PetscClassId id; 686 void (*func)(); 687 PetscInt numids; 688 const PetscInt *ids; 689 void *ctx; 690 691 ierr = DMGetBoundary(dm, b, &isEssential, NULL, &labelname, &field, NULL, NULL, &func, &numids, &ids, &ctx);CHKERRQ(ierr); 692 if (!isEssential) continue; 693 ierr = DMGetLabel(dm, labelname, &label);CHKERRQ(ierr); 694 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 695 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 696 if (id == PETSCFE_CLASSID) { 697 ierr = DMPlexInsertBoundaryValues_FEM_Internal(dm, time, field, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);CHKERRQ(ierr); 698 } else if (id == PETSCFV_CLASSID) { 699 if (!faceGeomFVM) continue; 700 ierr = DMPlexInsertBoundaryValues_FVM_Internal(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, 701 field, label, numids, ids, (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);CHKERRQ(ierr); 702 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 703 } 704 PetscFunctionReturn(0); 705 } 706 707 #undef __FUNCT__ 708 #define __FUNCT__ "DMComputeL2Diff_Plex" 709 PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 710 { 711 const PetscInt debug = 0; 712 PetscSection section; 713 PetscQuadrature quad; 714 Vec localX; 715 PetscScalar *funcVal, *interpolant; 716 PetscReal *coords, *v0, *J, *invJ, detJ; 717 PetscReal localDiff = 0.0; 718 const PetscReal *quadPoints, *quadWeights; 719 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 720 PetscErrorCode ierr; 721 722 PetscFunctionBegin; 723 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 724 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 725 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 726 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 727 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 728 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 729 for (field = 0; field < numFields; ++field) { 730 PetscObject obj; 731 PetscClassId id; 732 PetscInt Nc; 733 734 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 735 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 736 if (id == PETSCFE_CLASSID) { 737 PetscFE fe = (PetscFE) obj; 738 739 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 740 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 741 } else if (id == PETSCFV_CLASSID) { 742 PetscFV fv = (PetscFV) obj; 743 744 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 745 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 746 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 747 numComponents += Nc; 748 } 749 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 750 ierr = DMProjectFunctionLocal(dm, 0.0, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 751 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 752 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 753 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 754 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 755 for (c = cStart; c < cEnd; ++c) { 756 PetscScalar *x = NULL; 757 PetscReal elemDiff = 0.0; 758 759 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 760 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 761 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 762 763 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 764 PetscObject obj; 765 PetscClassId id; 766 void * const ctx = ctxs ? ctxs[field] : NULL; 767 PetscInt Nb, Nc, q, fc; 768 769 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 770 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 771 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 772 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 773 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 774 if (debug) { 775 char title[1024]; 776 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 777 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 778 } 779 for (q = 0; q < numQuadPoints; ++q) { 780 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 781 ierr = (*funcs[field])(dim, time, coords, Nc, funcVal, ctx);CHKERRQ(ierr); 782 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 783 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 784 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 785 for (fc = 0; fc < Nc; ++fc) { 786 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ);CHKERRQ(ierr);} 787 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 788 } 789 } 790 fieldOffset += Nb*Nc; 791 } 792 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 793 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 794 localDiff += elemDiff; 795 } 796 ierr = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 797 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 798 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 799 *diff = PetscSqrtReal(*diff); 800 PetscFunctionReturn(0); 801 } 802 803 #undef __FUNCT__ 804 #define __FUNCT__ "DMPlexComputeL2GradientDiff" 805 /*@C 806 DMPlexComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h. 807 808 Input Parameters: 809 + dm - The DM 810 , time - The time 811 . funcs - The gradient functions to evaluate for each field component 812 . ctxs - Optional array of contexts to pass to each function, or NULL. 813 . X - The coefficient vector u_h 814 - n - The vector to project along 815 816 Output Parameter: 817 . diff - The diff ||(grad u - grad u_h) . n||_2 818 819 Level: developer 820 821 .seealso: DMProjectFunction(), DMComputeL2Diff() 822 @*/ 823 PetscErrorCode DMPlexComputeL2GradientDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff) 824 { 825 const PetscInt debug = 0; 826 PetscSection section; 827 PetscQuadrature quad; 828 Vec localX; 829 PetscScalar *funcVal, *interpolantVec; 830 PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ; 831 PetscReal localDiff = 0.0; 832 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset, comp; 833 PetscErrorCode ierr; 834 835 PetscFunctionBegin; 836 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 837 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 838 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 839 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 840 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 841 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 842 for (field = 0; field < numFields; ++field) { 843 PetscFE fe; 844 PetscInt Nc; 845 846 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 847 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 848 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 849 numComponents += Nc; 850 } 851 /* ierr = DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 852 ierr = PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr); 853 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 854 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 855 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 856 for (c = cStart; c < cEnd; ++c) { 857 PetscScalar *x = NULL; 858 PetscReal elemDiff = 0.0; 859 860 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 861 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 862 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 863 864 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 865 PetscFE fe; 866 void * const ctx = ctxs ? ctxs[field] : NULL; 867 const PetscReal *quadPoints, *quadWeights; 868 PetscReal *basisDer; 869 PetscInt numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g; 870 871 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 872 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 873 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 874 ierr = PetscFEGetNumComponents(fe, &Ncomp);CHKERRQ(ierr); 875 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 876 if (debug) { 877 char title[1024]; 878 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 879 ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr); 880 } 881 for (q = 0; q < numQuadPoints; ++q) { 882 for (d = 0; d < dim; d++) { 883 coords[d] = v0[d]; 884 for (e = 0; e < dim; e++) { 885 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 886 } 887 } 888 ierr = (*funcs[field])(dim, time, coords, n, numFields, funcVal, ctx);CHKERRQ(ierr); 889 for (fc = 0; fc < Ncomp; ++fc) { 890 PetscScalar interpolant = 0.0; 891 892 for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0; 893 for (f = 0; f < Nb; ++f) { 894 const PetscInt fidx = f*Ncomp+fc; 895 896 for (d = 0; d < dim; ++d) { 897 realSpaceDer[d] = 0.0; 898 for (g = 0; g < dim; ++g) { 899 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g]; 900 } 901 interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d]; 902 } 903 } 904 for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d]; 905 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);} 906 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 907 } 908 } 909 comp += Ncomp; 910 fieldOffset += Nb*Ncomp; 911 } 912 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 913 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 914 localDiff += elemDiff; 915 } 916 ierr = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr); 917 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 918 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 919 *diff = PetscSqrtReal(*diff); 920 PetscFunctionReturn(0); 921 } 922 923 #undef __FUNCT__ 924 #define __FUNCT__ "DMPlexComputeL2FieldDiff" 925 /*@C 926 DMPlexComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components. 927 928 Input Parameters: 929 + dm - The DM 930 . time - The time 931 . funcs - The functions to evaluate for each field component 932 . ctxs - Optional array of contexts to pass to each function, or NULL. 933 - X - The coefficient vector u_h 934 935 Output Parameter: 936 . diff - The array of differences, ||u^f - u^f_h||_2 937 938 Level: developer 939 940 .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2GradientDiff() 941 @*/ 942 PetscErrorCode DMPlexComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[]) 943 { 944 const PetscInt debug = 0; 945 PetscSection section; 946 PetscQuadrature quad; 947 Vec localX; 948 PetscScalar *funcVal, *interpolant; 949 PetscReal *coords, *v0, *J, *invJ, detJ; 950 PetscReal *localDiff; 951 const PetscReal *quadPoints, *quadWeights; 952 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 953 PetscErrorCode ierr; 954 955 PetscFunctionBegin; 956 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 957 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 958 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 959 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 960 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 961 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 962 for (field = 0; field < numFields; ++field) { 963 PetscObject obj; 964 PetscClassId id; 965 PetscInt Nc; 966 967 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 968 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 969 if (id == PETSCFE_CLASSID) { 970 PetscFE fe = (PetscFE) obj; 971 972 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 973 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 974 } else if (id == PETSCFV_CLASSID) { 975 PetscFV fv = (PetscFV) obj; 976 977 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 978 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 979 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 980 numComponents += Nc; 981 } 982 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 983 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 984 ierr = PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 985 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 986 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 987 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 988 for (c = cStart; c < cEnd; ++c) { 989 PetscScalar *x = NULL; 990 991 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 992 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 993 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 994 995 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 996 PetscObject obj; 997 PetscClassId id; 998 void * const ctx = ctxs ? ctxs[field] : NULL; 999 PetscInt Nb, Nc, q, fc; 1000 1001 PetscReal elemDiff = 0.0; 1002 1003 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1004 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1005 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1006 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1007 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1008 if (debug) { 1009 char title[1024]; 1010 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 1011 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 1012 } 1013 for (q = 0; q < numQuadPoints; ++q) { 1014 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 1015 ierr = (*funcs[field])(dim, time, coords, numFields, funcVal, ctx);CHKERRQ(ierr); 1016 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1017 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1018 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1019 for (fc = 0; fc < Nc; ++fc) { 1020 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);} 1021 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 1022 } 1023 } 1024 fieldOffset += Nb*Nc; 1025 localDiff[field] += elemDiff; 1026 } 1027 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1028 } 1029 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1030 ierr = MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1031 for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]); 1032 ierr = PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 1033 PetscFunctionReturn(0); 1034 } 1035 1036 #undef __FUNCT__ 1037 #define __FUNCT__ "DMPlexComputeL2DiffVec" 1038 /*@C 1039 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. 1040 1041 Input Parameters: 1042 + dm - The DM 1043 . time - The time 1044 . funcs - The functions to evaluate for each field component 1045 . ctxs - Optional array of contexts to pass to each function, or NULL. 1046 - X - The coefficient vector u_h 1047 1048 Output Parameter: 1049 . D - A Vec which holds the difference ||u - u_h||_2 for each cell 1050 1051 Level: developer 1052 1053 .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMPlexComputeL2GradientDiff() 1054 @*/ 1055 PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D) 1056 { 1057 PetscSection section; 1058 PetscQuadrature quad; 1059 Vec localX; 1060 PetscScalar *funcVal, *interpolant; 1061 PetscReal *coords, *v0, *J, *invJ, detJ; 1062 const PetscReal *quadPoints, *quadWeights; 1063 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 1064 PetscErrorCode ierr; 1065 1066 PetscFunctionBegin; 1067 ierr = VecSet(D, 0.0);CHKERRQ(ierr); 1068 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1069 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1070 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1071 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1072 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1073 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1074 for (field = 0; field < numFields; ++field) { 1075 PetscObject obj; 1076 PetscClassId id; 1077 PetscInt Nc; 1078 1079 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1080 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1081 if (id == PETSCFE_CLASSID) { 1082 PetscFE fe = (PetscFE) obj; 1083 1084 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1085 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1086 } else if (id == PETSCFV_CLASSID) { 1087 PetscFV fv = (PetscFV) obj; 1088 1089 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1090 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1091 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1092 numComponents += Nc; 1093 } 1094 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 1095 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 1096 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1097 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1098 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1099 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1100 for (c = cStart; c < cEnd; ++c) { 1101 PetscScalar *x = NULL; 1102 PetscScalar elemDiff = 0.0; 1103 1104 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1105 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 1106 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1107 1108 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1109 PetscObject obj; 1110 PetscClassId id; 1111 void * const ctx = ctxs ? ctxs[field] : NULL; 1112 PetscInt Nb, Nc, q, fc; 1113 1114 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1115 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1116 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1117 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1118 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1119 for (q = 0; q < numQuadPoints; ++q) { 1120 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 1121 ierr = (*funcs[field])(dim, time, coords, Nc, funcVal, ctx);CHKERRQ(ierr); 1122 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1123 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1124 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1125 for (fc = 0; fc < Nc; ++fc) { 1126 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 1127 } 1128 } 1129 fieldOffset += Nb*Nc; 1130 } 1131 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1132 ierr = VecSetValuesSection(D, section, c, &elemDiff, ADD_VALUES);CHKERRQ(ierr); 1133 } 1134 ierr = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 1135 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1136 ierr = VecSqrtAbs(D);CHKERRQ(ierr); 1137 PetscFunctionReturn(0); 1138 } 1139 1140 #undef __FUNCT__ 1141 #define __FUNCT__ "DMPlexComputeIntegralFEM" 1142 /*@ 1143 DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user 1144 1145 Input Parameters: 1146 + dm - The mesh 1147 . X - Local input vector 1148 - user - The user context 1149 1150 Output Parameter: 1151 . integral - Local integral for each field 1152 1153 Level: developer 1154 1155 .seealso: DMPlexComputeResidualFEM() 1156 @*/ 1157 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user) 1158 { 1159 DM_Plex *mesh = (DM_Plex *) dm->data; 1160 DM dmAux; 1161 Vec localX, A; 1162 PetscDS prob, probAux = NULL; 1163 PetscSection section, sectionAux; 1164 PetscFECellGeom *cgeom; 1165 PetscScalar *u, *a = NULL; 1166 PetscReal *lintegral, *vol; 1167 PetscInt dim, Nf, f, numCells, cStart, cEnd, cEndInterior, c; 1168 PetscInt totDim, totDimAux; 1169 PetscErrorCode ierr; 1170 1171 PetscFunctionBegin; 1172 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1173 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1174 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1175 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1176 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1177 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1178 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1179 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1180 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1181 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1182 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1183 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1184 numCells = cEnd - cStart; 1185 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1186 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1187 if (dmAux) { 1188 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1189 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1190 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1191 } 1192 ierr = DMPlexInsertBoundaryValues(dm, localX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 1193 ierr = PetscMalloc4(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeom,numCells,&vol);CHKERRQ(ierr); 1194 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1195 for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;} 1196 for (c = cStart; c < cEnd; ++c) { 1197 PetscScalar *x = NULL; 1198 PetscInt i; 1199 1200 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeom[c].v0, cgeom[c].J, cgeom[c].invJ, &cgeom[c].detJ);CHKERRQ(ierr); 1201 ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol[c], NULL, NULL);CHKERRQ(ierr); 1202 if (cgeom[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeom[c].detJ, c); 1203 ierr = DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1204 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 1205 ierr = DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1206 if (dmAux) { 1207 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1208 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 1209 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1210 } 1211 } 1212 for (f = 0; f < Nf; ++f) { 1213 PetscObject obj; 1214 PetscClassId id; 1215 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 1216 1217 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1218 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1219 if (id == PETSCFE_CLASSID) { 1220 PetscFE fe = (PetscFE) obj; 1221 PetscQuadrature q; 1222 PetscInt Nq, Nb; 1223 1224 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1225 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1226 ierr = PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 1227 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1228 blockSize = Nb*Nq; 1229 batchSize = numBlocks * blockSize; 1230 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1231 numChunks = numCells / (numBatches*batchSize); 1232 Ne = numChunks*numBatches*batchSize; 1233 Nr = numCells % (numBatches*batchSize); 1234 offset = numCells - Nr; 1235 ierr = PetscFEIntegrate(fe, prob, f, Ne, cgeom, u, probAux, a, lintegral);CHKERRQ(ierr); 1236 ierr = PetscFEIntegrate(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);CHKERRQ(ierr); 1237 } else if (id == PETSCFV_CLASSID) { 1238 /* PetscFV fv = (PetscFV) obj; */ 1239 PetscInt foff; 1240 PetscPointFunc obj_func; 1241 PetscScalar lint; 1242 1243 ierr = PetscDSGetObjective(prob, f, &obj_func);CHKERRQ(ierr); 1244 ierr = PetscDSGetFieldOffset(prob, f, &foff);CHKERRQ(ierr); 1245 if (obj_func) { 1246 for (c = 0; c < numCells; ++c) { 1247 /* TODO: Need full pointwise interpolation and get centroid for x argument */ 1248 obj_func(dim, Nf, 0, NULL, NULL, &u[totDim*c+foff], NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0.0, NULL, &lint); 1249 lintegral[f] = PetscRealPart(lint)*vol[c]; 1250 } 1251 } 1252 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 1253 } 1254 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1255 if (mesh->printFEM) { 1256 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");CHKERRQ(ierr); 1257 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);CHKERRQ(ierr);} 1258 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 1259 } 1260 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1261 ierr = MPIU_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr); 1262 ierr = PetscFree4(lintegral,u,cgeom,vol);CHKERRQ(ierr); 1263 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1264 PetscFunctionReturn(0); 1265 } 1266 1267 #undef __FUNCT__ 1268 #define __FUNCT__ "DMPlexComputeInterpolatorNested" 1269 /*@ 1270 DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1271 1272 Input Parameters: 1273 + dmf - The fine mesh 1274 . dmc - The coarse mesh 1275 - user - The user context 1276 1277 Output Parameter: 1278 . In - The interpolation matrix 1279 1280 Level: developer 1281 1282 .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1283 @*/ 1284 PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user) 1285 { 1286 DM_Plex *mesh = (DM_Plex *) dmc->data; 1287 const char *name = "Interpolator"; 1288 PetscDS prob; 1289 PetscFE *feRef; 1290 PetscFV *fvRef; 1291 PetscSection fsection, fglobalSection; 1292 PetscSection csection, cglobalSection; 1293 PetscScalar *elemMat; 1294 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c; 1295 PetscInt cTotDim, rTotDim = 0; 1296 PetscErrorCode ierr; 1297 1298 PetscFunctionBegin; 1299 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1300 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1301 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1302 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1303 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1304 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1305 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1306 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1307 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1308 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1309 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr); 1310 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1311 for (f = 0; f < Nf; ++f) { 1312 PetscObject obj; 1313 PetscClassId id; 1314 PetscInt rNb = 0, Nc = 0; 1315 1316 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1317 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1318 if (id == PETSCFE_CLASSID) { 1319 PetscFE fe = (PetscFE) obj; 1320 1321 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1322 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1323 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1324 } else if (id == PETSCFV_CLASSID) { 1325 PetscFV fv = (PetscFV) obj; 1326 PetscDualSpace Q; 1327 1328 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1329 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1330 ierr = PetscDualSpaceGetDimension(Q, &rNb);CHKERRQ(ierr); 1331 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1332 } 1333 rTotDim += rNb*Nc; 1334 } 1335 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1336 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 1337 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 1338 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1339 PetscDualSpace Qref; 1340 PetscQuadrature f; 1341 const PetscReal *qpoints, *qweights; 1342 PetscReal *points; 1343 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1344 1345 /* Compose points from all dual basis functionals */ 1346 if (feRef[fieldI]) { 1347 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1348 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 1349 } else { 1350 ierr = PetscFVGetDualSpace(fvRef[fieldI], &Qref);CHKERRQ(ierr); 1351 ierr = PetscFVGetNumComponents(fvRef[fieldI], &Nc);CHKERRQ(ierr); 1352 } 1353 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1354 for (i = 0; i < fpdim; ++i) { 1355 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1356 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1357 npoints += Np; 1358 } 1359 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1360 for (i = 0, k = 0; i < fpdim; ++i) { 1361 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1362 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1363 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1364 } 1365 1366 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1367 PetscObject obj; 1368 PetscClassId id; 1369 PetscReal *B; 1370 PetscInt NcJ = 0, cpdim = 0, j; 1371 1372 ierr = PetscDSGetDiscretization(prob, fieldJ, &obj);CHKERRQ(ierr); 1373 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1374 if (id == PETSCFE_CLASSID) { 1375 PetscFE fe = (PetscFE) obj; 1376 1377 /* Evaluate basis at points */ 1378 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 1379 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1380 /* For now, fields only interpolate themselves */ 1381 if (fieldI == fieldJ) { 1382 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); 1383 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1384 for (i = 0, k = 0; i < fpdim; ++i) { 1385 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1386 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1387 for (p = 0; p < Np; ++p, ++k) { 1388 for (j = 0; j < cpdim; ++j) { 1389 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]; 1390 } 1391 } 1392 } 1393 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1394 } 1395 } else if (id == PETSCFV_CLASSID) { 1396 PetscFV fv = (PetscFV) obj; 1397 1398 /* Evaluate constant function at points */ 1399 ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr); 1400 cpdim = 1; 1401 /* For now, fields only interpolate themselves */ 1402 if (fieldI == fieldJ) { 1403 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); 1404 for (i = 0, k = 0; i < fpdim; ++i) { 1405 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1406 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1407 for (p = 0; p < Np; ++p, ++k) { 1408 for (j = 0; j < cpdim; ++j) { 1409 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p]; 1410 } 1411 } 1412 } 1413 } 1414 } 1415 offsetJ += cpdim*NcJ; 1416 } 1417 offsetI += fpdim*Nc; 1418 ierr = PetscFree(points);CHKERRQ(ierr); 1419 } 1420 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1421 /* Preallocate matrix */ 1422 { 1423 PetscHashJK ht; 1424 PetscLayout rLayout; 1425 PetscInt *dnz, *onz, *cellCIndices, *cellFIndices; 1426 PetscInt locRows, rStart, rEnd, cell, r; 1427 1428 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1429 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1430 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1431 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1432 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1433 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1434 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1435 ierr = PetscCalloc4(locRows,&dnz,locRows,&onz,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1436 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1437 for (cell = cStart; cell < cEnd; ++cell) { 1438 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1439 for (r = 0; r < rTotDim; ++r) { 1440 PetscHashJKKey key; 1441 PetscHashJKIter missing, iter; 1442 1443 key.j = cellFIndices[r]; 1444 if (key.j < 0) continue; 1445 for (c = 0; c < cTotDim; ++c) { 1446 key.k = cellCIndices[c]; 1447 if (key.k < 0) continue; 1448 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1449 if (missing) { 1450 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1451 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1452 else ++onz[key.j-rStart]; 1453 } 1454 } 1455 } 1456 } 1457 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1458 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1459 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1460 ierr = PetscFree4(dnz,onz,cellCIndices,cellFIndices);CHKERRQ(ierr); 1461 } 1462 /* Fill matrix */ 1463 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1464 for (c = cStart; c < cEnd; ++c) { 1465 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1466 } 1467 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1468 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1469 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1470 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1471 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1472 if (mesh->printFEM) { 1473 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1474 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1475 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1476 } 1477 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1478 PetscFunctionReturn(0); 1479 } 1480 1481 #undef __FUNCT__ 1482 #define __FUNCT__ "DMPlexComputeInterpolatorGeneral" 1483 /*@ 1484 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 1485 1486 Input Parameters: 1487 + dmf - The fine mesh 1488 . dmc - The coarse mesh 1489 - user - The user context 1490 1491 Output Parameter: 1492 . In - The interpolation matrix 1493 1494 Level: developer 1495 1496 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1497 @*/ 1498 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 1499 { 1500 PetscDS prob; 1501 PetscSection fsection, csection, globalFSection, globalCSection; 1502 PetscHashJK ht; 1503 PetscLayout rLayout; 1504 PetscInt *dnz, *onz; 1505 PetscInt locRows, rStart, rEnd; 1506 PetscReal *x, *v0, *J, *invJ, detJ; 1507 PetscReal *v0c, *Jc, *invJc, detJc; 1508 PetscScalar *elemMat; 1509 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1510 PetscErrorCode ierr; 1511 1512 PetscFunctionBegin; 1513 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1514 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1515 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1516 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1517 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1518 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1519 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1520 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1521 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1522 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1523 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1524 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1525 ierr = PetscMalloc1(totDim*totDim, &elemMat);CHKERRQ(ierr); 1526 1527 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1528 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1529 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1530 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1531 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1532 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1533 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1534 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1535 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1536 for (field = 0; field < Nf; ++field) { 1537 PetscObject obj; 1538 PetscClassId id; 1539 PetscDualSpace Q = NULL; 1540 PetscQuadrature f; 1541 const PetscReal *qpoints; 1542 PetscInt Nc, Np, fpdim, i, d; 1543 1544 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1545 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1546 if (id == PETSCFE_CLASSID) { 1547 PetscFE fe = (PetscFE) obj; 1548 1549 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1550 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1551 } else if (id == PETSCFV_CLASSID) { 1552 PetscFV fv = (PetscFV) obj; 1553 1554 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1555 Nc = 1; 1556 } 1557 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1558 /* For each fine grid cell */ 1559 for (cell = cStart; cell < cEnd; ++cell) { 1560 PetscInt *findices, *cindices; 1561 PetscInt numFIndices, numCIndices; 1562 1563 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1564 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1565 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1566 for (i = 0; i < fpdim; ++i) { 1567 Vec pointVec; 1568 PetscScalar *pV; 1569 IS coarseCellIS; 1570 const PetscInt *coarseCells; 1571 PetscInt numCoarseCells, q, r, c; 1572 1573 /* Get points from the dual basis functional quadrature */ 1574 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1575 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1576 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1577 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1578 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1579 for (q = 0; q < Np; ++q) { 1580 /* Transform point to real space */ 1581 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1582 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1583 } 1584 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1585 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1586 ierr = DMLocatePoints(dmc, pointVec, &coarseCellIS);CHKERRQ(ierr); 1587 ierr = ISViewFromOptions(coarseCellIS, NULL, "-interp_is_view");CHKERRQ(ierr); 1588 /* Update preallocation info */ 1589 ierr = ISGetLocalSize(coarseCellIS, &numCoarseCells);CHKERRQ(ierr); 1590 ierr = ISGetIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1591 for (r = 0; r < Nc; ++r) { 1592 PetscHashJKKey key; 1593 PetscHashJKIter missing, iter; 1594 1595 key.j = findices[i*Nc+r]; 1596 if (key.j < 0) continue; 1597 /* Get indices for coarse elements */ 1598 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1599 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1600 for (c = 0; c < numCIndices; ++c) { 1601 key.k = cindices[c]; 1602 if (key.k < 0) continue; 1603 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1604 if (missing) { 1605 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1606 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1607 else ++onz[key.j-rStart]; 1608 } 1609 } 1610 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1611 } 1612 } 1613 ierr = ISRestoreIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1614 ierr = ISDestroy(&coarseCellIS);CHKERRQ(ierr); 1615 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1616 } 1617 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1618 } 1619 } 1620 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1621 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1622 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1623 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1624 for (field = 0; field < Nf; ++field) { 1625 PetscObject obj; 1626 PetscClassId id; 1627 PetscDualSpace Q = NULL; 1628 PetscQuadrature f; 1629 const PetscReal *qpoints, *qweights; 1630 PetscInt Nc, Np, fpdim, i, d; 1631 1632 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1633 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1634 if (id == PETSCFE_CLASSID) { 1635 PetscFE fe = (PetscFE) obj; 1636 1637 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1638 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1639 } else if (id == PETSCFV_CLASSID) { 1640 PetscFV fv = (PetscFV) obj; 1641 1642 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1643 Nc = 1; 1644 } 1645 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1646 /* For each fine grid cell */ 1647 for (cell = cStart; cell < cEnd; ++cell) { 1648 PetscInt *findices, *cindices; 1649 PetscInt numFIndices, numCIndices; 1650 1651 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1652 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1653 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1654 for (i = 0; i < fpdim; ++i) { 1655 Vec pointVec; 1656 PetscScalar *pV; 1657 IS coarseCellIS; 1658 const PetscInt *coarseCells; 1659 PetscInt numCoarseCells, cpdim, q, c, j; 1660 1661 /* Get points from the dual basis functional quadrature */ 1662 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1663 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, &qweights);CHKERRQ(ierr); 1664 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1665 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1666 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1667 for (q = 0; q < Np; ++q) { 1668 /* Transform point to real space */ 1669 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1670 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1671 } 1672 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1673 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1674 ierr = DMLocatePoints(dmc, pointVec, &coarseCellIS);CHKERRQ(ierr); 1675 /* Update preallocation info */ 1676 ierr = ISGetLocalSize(coarseCellIS, &numCoarseCells);CHKERRQ(ierr); 1677 ierr = ISGetIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1678 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1679 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1680 PetscReal pVReal[3]; 1681 1682 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1683 /* Transform points from real space to coarse reference space */ 1684 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell], NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 1685 CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x); 1686 for (d = 0; d < dim; ++d) pV[ccell*dim+d] = pVReal[d]; 1687 1688 if (id == PETSCFE_CLASSID) { 1689 PetscFE fe = (PetscFE) obj; 1690 PetscReal *B; 1691 1692 /* Evaluate coarse basis on contained point */ 1693 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1694 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 1695 /* Get elemMat entries by multiplying by weight */ 1696 for (j = 0; j < cpdim; ++j) { 1697 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = B[j*Nc + c]*qweights[ccell]; 1698 } 1699 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1700 } else { 1701 cpdim = 1; 1702 for (j = 0; j < cpdim; ++j) { 1703 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = 1.0*qweights[ccell]; 1704 } 1705 } 1706 /* Update interpolator */ 1707 ierr = MatSetValues(In, Nc, &findices[i*Nc], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1708 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1709 } 1710 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1711 ierr = ISRestoreIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1712 ierr = ISDestroy(&coarseCellIS);CHKERRQ(ierr); 1713 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1714 } 1715 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1716 } 1717 } 1718 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 1719 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 1720 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1721 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1722 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1723 PetscFunctionReturn(0); 1724 } 1725 1726 #undef __FUNCT__ 1727 #define __FUNCT__ "DMPlexComputeInjectorFEM" 1728 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 1729 { 1730 PetscDS prob; 1731 PetscFE *feRef; 1732 PetscFV *fvRef; 1733 Vec fv, cv; 1734 IS fis, cis; 1735 PetscSection fsection, fglobalSection, csection, cglobalSection; 1736 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 1737 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, offsetC, offsetF, m; 1738 PetscErrorCode ierr; 1739 1740 PetscFunctionBegin; 1741 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1742 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1743 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1744 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1745 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1746 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1747 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1748 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1749 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1750 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1751 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1752 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1753 for (f = 0; f < Nf; ++f) { 1754 PetscObject obj; 1755 PetscClassId id; 1756 PetscInt fNb = 0, Nc = 0; 1757 1758 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1759 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1760 if (id == PETSCFE_CLASSID) { 1761 PetscFE fe = (PetscFE) obj; 1762 1763 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1764 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 1765 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1766 } else if (id == PETSCFV_CLASSID) { 1767 PetscFV fv = (PetscFV) obj; 1768 PetscDualSpace Q; 1769 1770 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1771 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1772 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 1773 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1774 } 1775 fTotDim += fNb*Nc; 1776 } 1777 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1778 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 1779 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 1780 PetscFE feC; 1781 PetscFV fvC; 1782 PetscDualSpace QF, QC; 1783 PetscInt NcF, NcC, fpdim, cpdim; 1784 1785 if (feRef[field]) { 1786 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 1787 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 1788 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 1789 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 1790 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1791 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 1792 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1793 } else { 1794 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 1795 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 1796 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 1797 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 1798 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1799 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 1800 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1801 } 1802 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); 1803 for (c = 0; c < cpdim; ++c) { 1804 PetscQuadrature cfunc; 1805 const PetscReal *cqpoints; 1806 PetscInt NpC; 1807 PetscBool found = PETSC_FALSE; 1808 1809 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 1810 ierr = PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 1811 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 1812 for (f = 0; f < fpdim; ++f) { 1813 PetscQuadrature ffunc; 1814 const PetscReal *fqpoints; 1815 PetscReal sum = 0.0; 1816 PetscInt NpF, comp; 1817 1818 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 1819 ierr = PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 1820 if (NpC != NpF) continue; 1821 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 1822 if (sum > 1.0e-9) continue; 1823 for (comp = 0; comp < NcC; ++comp) { 1824 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 1825 } 1826 found = PETSC_TRUE; 1827 break; 1828 } 1829 if (!found) { 1830 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 1831 if (fvRef[field]) { 1832 PetscInt comp; 1833 for (comp = 0; comp < NcC; ++comp) { 1834 cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp; 1835 } 1836 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 1837 } 1838 } 1839 offsetC += cpdim*NcC; 1840 offsetF += fpdim*NcF; 1841 } 1842 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 1843 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1844 1845 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 1846 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 1847 ierr = VecGetOwnershipRange(cv, &startC, NULL);CHKERRQ(ierr); 1848 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 1849 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 1850 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 1851 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 1852 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 1853 for (c = cStart; c < cEnd; ++c) { 1854 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 1855 for (d = 0; d < cTotDim; ++d) { 1856 if (cellCIndices[d] < 0) continue; 1857 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]]); 1858 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 1859 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 1860 } 1861 } 1862 ierr = PetscFree(cmap);CHKERRQ(ierr); 1863 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 1864 1865 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 1866 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 1867 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 1868 ierr = ISDestroy(&cis);CHKERRQ(ierr); 1869 ierr = ISDestroy(&fis);CHKERRQ(ierr); 1870 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 1871 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 1872 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1873 PetscFunctionReturn(0); 1874 } 1875