1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petscsf.h> 3 #include <petscsnes.h> 4 5 #include <petsc/private/hashsetij.h> 6 #include <petsc/private/petscfeimpl.h> 7 #include <petsc/private/petscfvimpl.h> 8 9 static PetscErrorCode DMPlexConvertPlex(DM dm, DM *plex, PetscBool copy) 10 { 11 PetscBool isPlex; 12 PetscErrorCode ierr; 13 14 PetscFunctionBegin; 15 ierr = PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);CHKERRQ(ierr); 16 if (isPlex) { 17 *plex = dm; 18 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 19 } else { 20 ierr = PetscObjectQuery((PetscObject) dm, "dm_plex", (PetscObject *) plex);CHKERRQ(ierr); 21 if (!*plex) { 22 ierr = DMConvert(dm,DMPLEX,plex);CHKERRQ(ierr); 23 ierr = PetscObjectCompose((PetscObject) dm, "dm_plex", (PetscObject) *plex);CHKERRQ(ierr); 24 if (copy) { 25 const char *comps[3] = {"A", "dmAux"}; 26 PetscObject obj; 27 PetscInt i; 28 29 { 30 /* Run the subdomain hook (this will copy the DMSNES/DMTS) */ 31 DMSubDomainHookLink link; 32 for (link = dm->subdomainhook; link; link = link->next) { 33 if (link->ddhook) {ierr = (*link->ddhook)(dm, *plex, link->ctx);CHKERRQ(ierr);} 34 } 35 } 36 for (i = 0; i < 3; i++) { 37 ierr = PetscObjectQuery((PetscObject) dm, comps[i], &obj);CHKERRQ(ierr); 38 ierr = PetscObjectCompose((PetscObject) *plex, comps[i], obj);CHKERRQ(ierr); 39 } 40 } 41 } else { 42 ierr = PetscObjectReference((PetscObject) *plex);CHKERRQ(ierr); 43 } 44 } 45 PetscFunctionReturn(0); 46 } 47 48 static PetscErrorCode PetscContainerUserDestroy_PetscFEGeom (void *ctx) 49 { 50 PetscFEGeom *geom = (PetscFEGeom *) ctx; 51 PetscErrorCode ierr; 52 53 PetscFunctionBegin; 54 ierr = PetscFEGeomDestroy(&geom);CHKERRQ(ierr); 55 PetscFunctionReturn(0); 56 } 57 58 static PetscErrorCode DMPlexGetFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom) 59 { 60 char composeStr[33] = {0}; 61 PetscObjectId id; 62 PetscContainer container; 63 PetscErrorCode ierr; 64 65 PetscFunctionBegin; 66 ierr = PetscObjectGetId((PetscObject)quad,&id);CHKERRQ(ierr); 67 ierr = PetscSNPrintf(composeStr, 32, "DMPlexGetFEGeom_%x\n", id);CHKERRQ(ierr); 68 ierr = PetscObjectQuery((PetscObject) pointIS, composeStr, (PetscObject *) &container);CHKERRQ(ierr); 69 if (container) { 70 ierr = PetscContainerGetPointer(container, (void **) geom);CHKERRQ(ierr); 71 } else { 72 ierr = DMFieldCreateFEGeom(coordField, pointIS, quad, faceData, geom);CHKERRQ(ierr); 73 ierr = PetscContainerCreate(PETSC_COMM_SELF,&container);CHKERRQ(ierr); 74 ierr = PetscContainerSetPointer(container, (void *) *geom);CHKERRQ(ierr); 75 ierr = PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom);CHKERRQ(ierr); 76 ierr = PetscObjectCompose((PetscObject) pointIS, composeStr, (PetscObject) container);CHKERRQ(ierr); 77 ierr = PetscContainerDestroy(&container);CHKERRQ(ierr); 78 } 79 PetscFunctionReturn(0); 80 } 81 82 static PetscErrorCode DMPlexRestoreFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom) 83 { 84 PetscFunctionBegin; 85 *geom = NULL; 86 PetscFunctionReturn(0); 87 } 88 89 /*@ 90 DMPlexGetScale - Get the scale for the specified fundamental unit 91 92 Not collective 93 94 Input Arguments: 95 + dm - the DM 96 - unit - The SI unit 97 98 Output Argument: 99 . scale - The value used to scale all quantities with this unit 100 101 Level: advanced 102 103 .seealso: DMPlexSetScale(), PetscUnit 104 @*/ 105 PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale) 106 { 107 DM_Plex *mesh = (DM_Plex*) dm->data; 108 109 PetscFunctionBegin; 110 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 111 PetscValidPointer(scale, 3); 112 *scale = mesh->scale[unit]; 113 PetscFunctionReturn(0); 114 } 115 116 /*@ 117 DMPlexSetScale - Set the scale for the specified fundamental unit 118 119 Not collective 120 121 Input Arguments: 122 + dm - the DM 123 . unit - The SI unit 124 - scale - The value used to scale all quantities with this unit 125 126 Level: advanced 127 128 .seealso: DMPlexGetScale(), PetscUnit 129 @*/ 130 PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale) 131 { 132 DM_Plex *mesh = (DM_Plex*) dm->data; 133 134 PetscFunctionBegin; 135 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 136 mesh->scale[unit] = scale; 137 PetscFunctionReturn(0); 138 } 139 140 static PetscErrorCode DMPlexProjectRigidBody_Private(PetscInt dim, PetscReal t, const PetscReal X[], PetscInt Nf, PetscScalar *mode, void *ctx) 141 { 142 const PetscInt eps[3][3][3] = {{{0, 0, 0}, {0, 0, 1}, {0, -1, 0}}, {{0, 0, -1}, {0, 0, 0}, {1, 0, 0}}, {{0, 1, 0}, {-1, 0, 0}, {0, 0, 0}}}; 143 PetscInt *ctxInt = (PetscInt *) ctx; 144 PetscInt dim2 = ctxInt[0]; 145 PetscInt d = ctxInt[1]; 146 PetscInt i, j, k = dim > 2 ? d - dim : d; 147 148 PetscFunctionBegin; 149 if (dim != dim2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input dimension %d does not match context dimension %d", dim, dim2); 150 for (i = 0; i < dim; i++) mode[i] = 0.; 151 if (d < dim) { 152 mode[d] = 1.; /* Translation along axis d */ 153 } else { 154 for (i = 0; i < dim; i++) { 155 for (j = 0; j < dim; j++) { 156 mode[j] += eps[i][j][k]*X[i]; /* Rotation about axis d */ 157 } 158 } 159 } 160 PetscFunctionReturn(0); 161 } 162 163 /*@ 164 DMPlexCreateRigidBody - For the default global section, create rigid body modes by function space interpolation 165 166 Collective on DM 167 168 Input Arguments: 169 . dm - the DM 170 171 Output Argument: 172 . sp - the null space 173 174 Note: This is necessary to provide a suitable coarse space for algebraic multigrid 175 176 Level: advanced 177 178 .seealso: MatNullSpaceCreate(), PCGAMG 179 @*/ 180 PetscErrorCode DMPlexCreateRigidBody(DM dm, MatNullSpace *sp) 181 { 182 MPI_Comm comm; 183 Vec mode[6]; 184 PetscSection section, globalSection; 185 PetscInt dim, dimEmbed, n, m, mmin, d, i, j; 186 PetscErrorCode ierr; 187 188 PetscFunctionBegin; 189 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 190 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 191 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 192 if (dim == 1) { 193 ierr = MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);CHKERRQ(ierr); 194 PetscFunctionReturn(0); 195 } 196 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 197 ierr = DMGetGlobalSection(dm, &globalSection);CHKERRQ(ierr); 198 ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr); 199 m = (dim*(dim+1))/2; 200 ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr); 201 ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr); 202 ierr = VecSetUp(mode[0]);CHKERRQ(ierr); 203 ierr = VecGetSize(mode[0], &n);CHKERRQ(ierr); 204 mmin = PetscMin(m, n); 205 for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);} 206 for (d = 0; d < m; d++) { 207 PetscInt ctx[2]; 208 PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private; 209 void *voidctx = (void *) (&ctx[0]); 210 211 ctx[0] = dimEmbed; 212 ctx[1] = d; 213 ierr = DMProjectFunction(dm, 0.0, &func, &voidctx, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 214 } 215 for (i = 0; i < PetscMin(dim, mmin); ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);} 216 /* Orthonormalize system */ 217 for (i = dim; i < mmin; ++i) { 218 PetscScalar dots[6]; 219 220 ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr); 221 for (j = 0; j < i; ++j) dots[j] *= -1.0; 222 ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr); 223 ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr); 224 } 225 ierr = MatNullSpaceCreate(comm, PETSC_FALSE, mmin, mode, sp);CHKERRQ(ierr); 226 for (i = 0; i < m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);} 227 PetscFunctionReturn(0); 228 } 229 230 /*@ 231 DMPlexCreateRigidBodies - For the default global section, create rigid body modes by function space interpolation 232 233 Collective on DM 234 235 Input Arguments: 236 + dm - the DM 237 . nb - The number of bodies 238 . label - The DMLabel marking each domain 239 . nids - The number of ids per body 240 - ids - An array of the label ids in sequence for each domain 241 242 Output Argument: 243 . sp - the null space 244 245 Note: This is necessary to provide a suitable coarse space for algebraic multigrid 246 247 Level: advanced 248 249 .seealso: MatNullSpaceCreate() 250 @*/ 251 PetscErrorCode DMPlexCreateRigidBodies(DM dm, PetscInt nb, DMLabel label, const PetscInt nids[], const PetscInt ids[], MatNullSpace *sp) 252 { 253 MPI_Comm comm; 254 PetscSection section, globalSection; 255 Vec *mode; 256 PetscScalar *dots; 257 PetscInt dim, dimEmbed, n, m, b, d, i, j, off; 258 PetscErrorCode ierr; 259 260 PetscFunctionBegin; 261 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 262 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 263 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 264 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 265 ierr = DMGetGlobalSection(dm, &globalSection);CHKERRQ(ierr); 266 ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr); 267 m = nb * (dim*(dim+1))/2; 268 ierr = PetscMalloc2(m, &mode, m, &dots);CHKERRQ(ierr); 269 ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr); 270 ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr); 271 ierr = VecSetUp(mode[0]);CHKERRQ(ierr); 272 for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);} 273 for (b = 0, off = 0; b < nb; ++b) { 274 for (d = 0; d < m/nb; ++d) { 275 PetscInt ctx[2]; 276 PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private; 277 void *voidctx = (void *) (&ctx[0]); 278 279 ctx[0] = dimEmbed; 280 ctx[1] = d; 281 ierr = DMProjectFunctionLabel(dm, 0.0, label, nids[b], &ids[off], 0, NULL, &func, &voidctx, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 282 off += nids[b]; 283 } 284 } 285 for (i = 0; i < dim; ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);} 286 /* Orthonormalize system */ 287 for (i = 0; i < m; ++i) { 288 ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr); 289 for (j = 0; j < i; ++j) dots[j] *= -1.0; 290 ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr); 291 ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr); 292 } 293 ierr = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);CHKERRQ(ierr); 294 for (i = 0; i< m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);} 295 ierr = PetscFree2(mode, dots);CHKERRQ(ierr); 296 PetscFunctionReturn(0); 297 } 298 299 /*@ 300 DMPlexSetMaxProjectionHeight - In DMPlexProjectXXXLocal() functions, the projected values of a basis function's dofs 301 are computed by associating the basis function with one of the mesh points in its transitively-closed support, and 302 evaluating the dual space basis of that point. A basis function is associated with the point in its 303 transitively-closed support whose mesh height is highest (w.r.t. DAG height), but not greater than the maximum 304 projection height, which is set with this function. By default, the maximum projection height is zero, which means 305 that only mesh cells are used to project basis functions. A height of one, for example, evaluates a cell-interior 306 basis functions using its cells dual space basis, but all other basis functions with the dual space basis of a face. 307 308 Input Parameters: 309 + dm - the DMPlex object 310 - height - the maximum projection height >= 0 311 312 Level: advanced 313 314 .seealso: DMPlexGetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal() 315 @*/ 316 PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height) 317 { 318 DM_Plex *plex = (DM_Plex *) dm->data; 319 320 PetscFunctionBegin; 321 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 322 plex->maxProjectionHeight = height; 323 PetscFunctionReturn(0); 324 } 325 326 /*@ 327 DMPlexGetMaxProjectionHeight - Get the maximum height (w.r.t. DAG) of mesh points used to evaluate dual bases in 328 DMPlexProjectXXXLocal() functions. 329 330 Input Parameters: 331 . dm - the DMPlex object 332 333 Output Parameters: 334 . height - the maximum projection height 335 336 Level: intermediate 337 338 .seealso: DMPlexSetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal() 339 @*/ 340 PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height) 341 { 342 DM_Plex *plex = (DM_Plex *) dm->data; 343 344 PetscFunctionBegin; 345 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 346 *height = plex->maxProjectionHeight; 347 PetscFunctionReturn(0); 348 } 349 350 /*@C 351 DMPlexInsertBoundaryValuesEssential - Insert boundary values into a local vector 352 353 Input Parameters: 354 + dm - The DM, with a PetscDS that matches the problem being constrained 355 . time - The time 356 . field - The field to constrain 357 . Nc - The number of constrained field components, or 0 for all components 358 . comps - An array of constrained component numbers, or NULL for all components 359 . label - The DMLabel defining constrained points 360 . numids - The number of DMLabel ids for constrained points 361 . ids - An array of ids for constrained points 362 . func - A pointwise function giving boundary values 363 - ctx - An optional user context for bcFunc 364 365 Output Parameter: 366 . locX - A local vector to receives the boundary values 367 368 Level: developer 369 370 .seealso: DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary() 371 @*/ 372 PetscErrorCode DMPlexInsertBoundaryValuesEssential(DM dm, PetscReal time, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void *ctx, Vec locX) 373 { 374 PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx); 375 void **ctxs; 376 PetscInt numFields; 377 PetscErrorCode ierr; 378 379 PetscFunctionBegin; 380 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 381 ierr = PetscCalloc2(numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 382 funcs[field] = func; 383 ctxs[field] = ctx; 384 ierr = DMProjectFunctionLabelLocal(dm, time, label, numids, ids, Nc, comps, funcs, ctxs, INSERT_BC_VALUES, locX);CHKERRQ(ierr); 385 ierr = PetscFree2(funcs,ctxs);CHKERRQ(ierr); 386 PetscFunctionReturn(0); 387 } 388 389 /*@C 390 DMPlexInsertBoundaryValuesEssentialField - Insert boundary values into a local vector 391 392 Input Parameters: 393 + dm - The DM, with a PetscDS that matches the problem being constrained 394 . time - The time 395 . locU - A local vector with the input solution values 396 . field - The field to constrain 397 . Nc - The number of constrained field components, or 0 for all components 398 . comps - An array of constrained component numbers, or NULL for all components 399 . label - The DMLabel defining constrained points 400 . numids - The number of DMLabel ids for constrained points 401 . ids - An array of ids for constrained points 402 . func - A pointwise function giving boundary values 403 - ctx - An optional user context for bcFunc 404 405 Output Parameter: 406 . locX - A local vector to receives the boundary values 407 408 Level: developer 409 410 .seealso: DMPlexInsertBoundaryValuesEssential(), DMAddBoundary() 411 @*/ 412 PetscErrorCode DMPlexInsertBoundaryValuesEssentialField(DM dm, PetscReal time, Vec locU, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[], 413 void (*func)(PetscInt, PetscInt, PetscInt, 414 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 415 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 416 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], 417 PetscScalar[]), 418 void *ctx, Vec locX) 419 { 420 void (**funcs)(PetscInt, PetscInt, PetscInt, 421 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 422 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 423 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]); 424 void **ctxs; 425 PetscInt numFields; 426 PetscErrorCode ierr; 427 428 PetscFunctionBegin; 429 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 430 ierr = PetscCalloc2(numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 431 funcs[field] = func; 432 ctxs[field] = ctx; 433 ierr = DMProjectFieldLabelLocal(dm, time, label, numids, ids, Nc, comps, locU, funcs, INSERT_BC_VALUES, locX);CHKERRQ(ierr); 434 ierr = PetscFree2(funcs,ctxs);CHKERRQ(ierr); 435 PetscFunctionReturn(0); 436 } 437 438 /*@C 439 DMPlexInsertBoundaryValuesRiemann - Insert boundary values into a local vector 440 441 Input Parameters: 442 + dm - The DM, with a PetscDS that matches the problem being constrained 443 . time - The time 444 . faceGeometry - A vector with the FVM face geometry information 445 . cellGeometry - A vector with the FVM cell geometry information 446 . Grad - A vector with the FVM cell gradient information 447 . field - The field to constrain 448 . Nc - The number of constrained field components, or 0 for all components 449 . comps - An array of constrained component numbers, or NULL for all components 450 . label - The DMLabel defining constrained points 451 . numids - The number of DMLabel ids for constrained points 452 . ids - An array of ids for constrained points 453 . func - A pointwise function giving boundary values 454 - ctx - An optional user context for bcFunc 455 456 Output Parameter: 457 . locX - A local vector to receives the boundary values 458 459 Note: This implementation currently ignores the numcomps/comps argument from DMAddBoundary() 460 461 Level: developer 462 463 .seealso: DMPlexInsertBoundaryValuesEssential(), DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary() 464 @*/ 465 PetscErrorCode DMPlexInsertBoundaryValuesRiemann(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[], 466 PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX) 467 { 468 PetscDS prob; 469 PetscSF sf; 470 DM dmFace, dmCell, dmGrad; 471 const PetscScalar *facegeom, *cellgeom = NULL, *grad; 472 const PetscInt *leaves; 473 PetscScalar *x, *fx; 474 PetscInt dim, nleaves, loc, fStart, fEnd, pdim, i; 475 PetscErrorCode ierr, ierru = 0; 476 477 PetscFunctionBegin; 478 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 479 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 480 nleaves = PetscMax(0, nleaves); 481 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 482 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 483 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 484 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 485 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 486 if (cellGeometry) { 487 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 488 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 489 } 490 if (Grad) { 491 PetscFV fv; 492 493 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);CHKERRQ(ierr); 494 ierr = VecGetDM(Grad, &dmGrad);CHKERRQ(ierr); 495 ierr = VecGetArrayRead(Grad, &grad);CHKERRQ(ierr); 496 ierr = PetscFVGetNumComponents(fv, &pdim);CHKERRQ(ierr); 497 ierr = DMGetWorkArray(dm, pdim, MPIU_SCALAR, &fx);CHKERRQ(ierr); 498 } 499 ierr = VecGetArray(locX, &x);CHKERRQ(ierr); 500 for (i = 0; i < numids; ++i) { 501 IS faceIS; 502 const PetscInt *faces; 503 PetscInt numFaces, f; 504 505 ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr); 506 if (!faceIS) continue; /* No points with that id on this process */ 507 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 508 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 509 for (f = 0; f < numFaces; ++f) { 510 const PetscInt face = faces[f], *cells; 511 PetscFVFaceGeom *fg; 512 513 if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */ 514 ierr = PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);CHKERRQ(ierr); 515 if (loc >= 0) continue; 516 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 517 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 518 if (Grad) { 519 PetscFVCellGeom *cg; 520 PetscScalar *cx, *cgrad; 521 PetscScalar *xG; 522 PetscReal dx[3]; 523 PetscInt d; 524 525 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);CHKERRQ(ierr); 526 ierr = DMPlexPointLocalRead(dm, cells[0], x, &cx);CHKERRQ(ierr); 527 ierr = DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);CHKERRQ(ierr); 528 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 529 DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx); 530 for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx); 531 ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx); 532 if (ierru) { 533 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 534 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 535 goto cleanup; 536 } 537 } else { 538 PetscScalar *xI; 539 PetscScalar *xG; 540 541 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr); 542 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 543 ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx); 544 if (ierru) { 545 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 546 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 547 goto cleanup; 548 } 549 } 550 } 551 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 552 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 553 } 554 cleanup: 555 ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr); 556 if (Grad) { 557 ierr = DMRestoreWorkArray(dm, pdim, MPIU_SCALAR, &fx);CHKERRQ(ierr); 558 ierr = VecRestoreArrayRead(Grad, &grad);CHKERRQ(ierr); 559 } 560 if (cellGeometry) {ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr);} 561 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 562 CHKERRQ(ierru); 563 PetscFunctionReturn(0); 564 } 565 566 PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM) 567 { 568 PetscDS prob; 569 PetscInt numBd, b; 570 PetscErrorCode ierr; 571 572 PetscFunctionBegin; 573 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 574 ierr = PetscDSGetNumBoundary(prob, &numBd);CHKERRQ(ierr); 575 for (b = 0; b < numBd; ++b) { 576 DMBoundaryConditionType type; 577 const char *labelname; 578 DMLabel label; 579 PetscInt field, Nc; 580 const PetscInt *comps; 581 PetscObject obj; 582 PetscClassId id; 583 void (*func)(void); 584 PetscInt numids; 585 const PetscInt *ids; 586 void *ctx; 587 588 ierr = DMGetBoundary(dm, b, &type, NULL, &labelname, &field, &Nc, &comps, &func, &numids, &ids, &ctx);CHKERRQ(ierr); 589 if (insertEssential != (type & DM_BC_ESSENTIAL)) continue; 590 ierr = DMGetLabel(dm, labelname, &label);CHKERRQ(ierr); 591 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 592 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 593 if (id == PETSCFE_CLASSID) { 594 switch (type) { 595 /* for FEM, there is no insertion to be done for non-essential boundary conditions */ 596 case DM_BC_ESSENTIAL: 597 ierr = DMPlexLabelAddCells(dm,label);CHKERRQ(ierr); 598 ierr = DMPlexInsertBoundaryValuesEssential(dm, time, field, Nc, comps, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);CHKERRQ(ierr); 599 ierr = DMPlexLabelClearCells(dm,label);CHKERRQ(ierr); 600 break; 601 case DM_BC_ESSENTIAL_FIELD: 602 ierr = DMPlexLabelAddCells(dm,label);CHKERRQ(ierr); 603 ierr = DMPlexInsertBoundaryValuesEssentialField(dm, time, locX, field, Nc, comps, label, numids, ids, 604 (void (*)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 605 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 606 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) func, ctx, locX);CHKERRQ(ierr); 607 ierr = DMPlexLabelClearCells(dm,label);CHKERRQ(ierr); 608 break; 609 default: break; 610 } 611 } else if (id == PETSCFV_CLASSID) { 612 if (!faceGeomFVM) continue; 613 ierr = DMPlexInsertBoundaryValuesRiemann(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, field, Nc, comps, label, numids, ids, 614 (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);CHKERRQ(ierr); 615 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 616 } 617 PetscFunctionReturn(0); 618 } 619 620 /*@ 621 DMPlexInsertBoundaryValues - Puts coefficients which represent boundary values into the local solution vector 622 623 Input Parameters: 624 + dm - The DM 625 . insertEssential - Should I insert essential (e.g. Dirichlet) or inessential (e.g. Neumann) boundary conditions 626 . time - The time 627 . faceGeomFVM - Face geometry data for FV discretizations 628 . cellGeomFVM - Cell geometry data for FV discretizations 629 - gradFVM - Gradient reconstruction data for FV discretizations 630 631 Output Parameters: 632 . locX - Solution updated with boundary values 633 634 Level: developer 635 636 .seealso: DMProjectFunctionLabelLocal() 637 @*/ 638 PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM) 639 { 640 PetscErrorCode ierr; 641 642 PetscFunctionBegin; 643 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 644 PetscValidHeaderSpecific(locX, VEC_CLASSID, 2); 645 if (faceGeomFVM) {PetscValidHeaderSpecific(faceGeomFVM, VEC_CLASSID, 4);} 646 if (cellGeomFVM) {PetscValidHeaderSpecific(cellGeomFVM, VEC_CLASSID, 5);} 647 if (gradFVM) {PetscValidHeaderSpecific(gradFVM, VEC_CLASSID, 6);} 648 ierr = PetscTryMethod(dm,"DMPlexInsertBoundaryValues_C",(DM,PetscBool,Vec,PetscReal,Vec,Vec,Vec),(dm,insertEssential,locX,time,faceGeomFVM,cellGeomFVM,gradFVM));CHKERRQ(ierr); 649 PetscFunctionReturn(0); 650 } 651 652 PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 653 { 654 Vec localX; 655 PetscErrorCode ierr; 656 657 PetscFunctionBegin; 658 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 659 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, time, NULL, NULL, NULL);CHKERRQ(ierr); 660 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 661 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 662 ierr = DMPlexComputeL2DiffLocal(dm, time, funcs, ctxs, localX, diff);CHKERRQ(ierr); 663 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 664 PetscFunctionReturn(0); 665 } 666 667 /*@C 668 DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h. 669 670 Input Parameters: 671 + dm - The DM 672 . time - The time 673 . funcs - The functions to evaluate for each field component 674 . ctxs - Optional array of contexts to pass to each function, or NULL. 675 - localX - The coefficient vector u_h, a local vector 676 677 Output Parameter: 678 . diff - The diff ||u - u_h||_2 679 680 Level: developer 681 682 .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff() 683 @*/ 684 PetscErrorCode DMPlexComputeL2DiffLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec localX, PetscReal *diff) 685 { 686 const PetscInt debug = ((DM_Plex*)dm->data)->printL2; 687 PetscSection section; 688 PetscQuadrature quad; 689 PetscScalar *funcVal, *interpolant; 690 PetscReal *coords, *detJ, *J; 691 PetscReal localDiff = 0.0; 692 const PetscReal *quadWeights; 693 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cellHeight, cStart, cEnd, cEndInterior, c, field, fieldOffset; 694 PetscErrorCode ierr; 695 696 PetscFunctionBegin; 697 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 698 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 699 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 700 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 701 for (field = 0; field < numFields; ++field) { 702 PetscObject obj; 703 PetscClassId id; 704 PetscInt Nc; 705 706 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 707 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 708 if (id == PETSCFE_CLASSID) { 709 PetscFE fe = (PetscFE) obj; 710 711 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 712 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 713 } else if (id == PETSCFV_CLASSID) { 714 PetscFV fv = (PetscFV) obj; 715 716 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 717 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 718 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 719 numComponents += Nc; 720 } 721 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);CHKERRQ(ierr); 722 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 723 ierr = PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);CHKERRQ(ierr); 724 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr); 725 ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr); 726 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 727 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 728 for (c = cStart; c < cEnd; ++c) { 729 PetscScalar *x = NULL; 730 PetscReal elemDiff = 0.0; 731 PetscInt qc = 0; 732 733 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);CHKERRQ(ierr); 734 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 735 736 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 737 PetscObject obj; 738 PetscClassId id; 739 void * const ctx = ctxs ? ctxs[field] : NULL; 740 PetscInt Nb, Nc, q, fc; 741 742 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 743 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 744 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 745 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 746 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 747 if (debug) { 748 char title[1024]; 749 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 750 ierr = DMPrintCellVector(c, title, Nb, &x[fieldOffset]);CHKERRQ(ierr); 751 } 752 for (q = 0; q < Nq; ++q) { 753 if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, point %D", detJ[q], c, q); 754 ierr = (*funcs[field])(coordDim, time, &coords[coordDim * q], Nc, funcVal, ctx); 755 if (ierr) { 756 PetscErrorCode ierr2; 757 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 758 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 759 ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2); 760 CHKERRQ(ierr); 761 } 762 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 763 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 764 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 765 for (fc = 0; fc < Nc; ++fc) { 766 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 767 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);CHKERRQ(ierr);} 768 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]; 769 } 770 } 771 fieldOffset += Nb; 772 qc += Nc; 773 } 774 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 775 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 776 localDiff += elemDiff; 777 } 778 ierr = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr); 779 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 780 *diff = PetscSqrtReal(*diff); 781 PetscFunctionReturn(0); 782 } 783 784 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) 785 { 786 const PetscInt debug = ((DM_Plex*)dm->data)->printL2; 787 PetscSection section; 788 PetscQuadrature quad; 789 Vec localX; 790 PetscScalar *funcVal, *interpolant; 791 const PetscReal *quadPoints, *quadWeights; 792 PetscReal *coords, *realSpaceDer, *J, *invJ, *detJ; 793 PetscReal localDiff = 0.0; 794 PetscInt dim, coordDim, qNc = 0, Nq = 0, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset; 795 PetscErrorCode ierr; 796 797 PetscFunctionBegin; 798 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 799 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 800 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 801 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 802 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 803 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 804 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 805 for (field = 0; field < numFields; ++field) { 806 PetscFE fe; 807 PetscInt Nc; 808 809 ierr = DMGetField(dm, field, NULL, (PetscObject *) &fe);CHKERRQ(ierr); 810 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 811 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 812 numComponents += Nc; 813 } 814 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 815 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 816 /* ierr = DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 817 ierr = PetscMalloc7(numComponents,&funcVal,coordDim*Nq,&coords,coordDim*Nq,&realSpaceDer,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ,numComponents,&interpolant,Nq,&detJ);CHKERRQ(ierr); 818 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 819 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 820 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 821 for (c = cStart; c < cEnd; ++c) { 822 PetscScalar *x = NULL; 823 PetscReal elemDiff = 0.0; 824 PetscInt qc = 0; 825 826 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, invJ, detJ);CHKERRQ(ierr); 827 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 828 829 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 830 PetscFE fe; 831 void * const ctx = ctxs ? ctxs[field] : NULL; 832 PetscReal *basisDer; 833 PetscInt Nb, Nc, q, fc; 834 835 ierr = DMGetField(dm, field, NULL, (PetscObject *) &fe);CHKERRQ(ierr); 836 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 837 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 838 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 839 if (debug) { 840 char title[1024]; 841 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 842 ierr = DMPrintCellVector(c, title, Nb, &x[fieldOffset]);CHKERRQ(ierr); 843 } 844 for (q = 0; q < Nq; ++q) { 845 if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", detJ[q], c, q); 846 ierr = (*funcs[field])(coordDim, time, &coords[q*coordDim], n, numFields, funcVal, ctx); 847 if (ierr) { 848 PetscErrorCode ierr2; 849 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 850 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 851 ierr2 = PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);CHKERRQ(ierr2); 852 CHKERRQ(ierr); 853 } 854 ierr = PetscFEInterpolateGradient_Static(fe, &x[fieldOffset], coordDim, invJ, n, q, interpolant);CHKERRQ(ierr); 855 for (fc = 0; fc < Nc; ++fc) { 856 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 857 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);CHKERRQ(ierr);} 858 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]; 859 } 860 } 861 fieldOffset += Nb; 862 qc += Nc; 863 } 864 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 865 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 866 localDiff += elemDiff; 867 } 868 ierr = PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);CHKERRQ(ierr); 869 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 870 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 871 *diff = PetscSqrtReal(*diff); 872 PetscFunctionReturn(0); 873 } 874 875 PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 876 { 877 const PetscInt debug = ((DM_Plex*)dm->data)->printL2; 878 PetscSection section; 879 PetscQuadrature quad; 880 Vec localX; 881 PetscScalar *funcVal, *interpolant; 882 PetscReal *coords, *detJ, *J; 883 PetscReal *localDiff; 884 const PetscReal *quadPoints, *quadWeights; 885 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset; 886 PetscErrorCode ierr; 887 888 PetscFunctionBegin; 889 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 890 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 891 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 892 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 893 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 894 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 895 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 896 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 897 for (field = 0; field < numFields; ++field) { 898 PetscObject obj; 899 PetscClassId id; 900 PetscInt Nc; 901 902 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 903 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 904 if (id == PETSCFE_CLASSID) { 905 PetscFE fe = (PetscFE) obj; 906 907 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 908 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 909 } else if (id == PETSCFV_CLASSID) { 910 PetscFV fv = (PetscFV) obj; 911 912 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 913 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 914 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 915 numComponents += Nc; 916 } 917 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 918 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 919 ierr = PetscCalloc6(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);CHKERRQ(ierr); 920 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 921 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 922 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 923 for (c = cStart; c < cEnd; ++c) { 924 PetscScalar *x = NULL; 925 PetscInt qc = 0; 926 927 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);CHKERRQ(ierr); 928 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 929 930 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 931 PetscObject obj; 932 PetscClassId id; 933 void * const ctx = ctxs ? ctxs[field] : NULL; 934 PetscInt Nb, Nc, q, fc; 935 936 PetscReal elemDiff = 0.0; 937 938 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 939 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 940 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 941 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 942 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 943 if (debug) { 944 char title[1024]; 945 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 946 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 947 } 948 for (q = 0; q < Nq; ++q) { 949 if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature point %D", detJ, c, q); 950 ierr = (*funcs[field])(coordDim, time, &coords[coordDim*q], numFields, funcVal, ctx); 951 if (ierr) { 952 PetscErrorCode ierr2; 953 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 954 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 955 ierr2 = PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2); 956 CHKERRQ(ierr); 957 } 958 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 959 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 960 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 961 for (fc = 0; fc < Nc; ++fc) { 962 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 963 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d field %d point %g %g %g diff %g\n", c, field, coordDim > 0 ? coords[0] : 0., coordDim > 1 ? coords[1] : 0., coordDim > 2 ? coords[2] : 0., PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);CHKERRQ(ierr);} 964 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]; 965 } 966 } 967 fieldOffset += Nb; 968 qc += Nc; 969 localDiff[field] += elemDiff; 970 } 971 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 972 } 973 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 974 ierr = MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 975 for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]); 976 ierr = PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr); 977 PetscFunctionReturn(0); 978 } 979 980 /*@C 981 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. 982 983 Input Parameters: 984 + dm - The DM 985 . time - The time 986 . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation 987 . ctxs - Optional array of contexts to pass to each function, or NULL. 988 - X - The coefficient vector u_h 989 990 Output Parameter: 991 . D - A Vec which holds the difference ||u - u_h||_2 for each cell 992 993 Level: developer 994 995 .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff() 996 @*/ 997 PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D) 998 { 999 PetscSection section; 1000 PetscQuadrature quad; 1001 Vec localX; 1002 PetscScalar *funcVal, *interpolant; 1003 PetscReal *coords, *detJ, *J; 1004 const PetscReal *quadPoints, *quadWeights; 1005 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset; 1006 PetscErrorCode ierr; 1007 1008 PetscFunctionBegin; 1009 ierr = VecSet(D, 0.0);CHKERRQ(ierr); 1010 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1011 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 1012 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 1013 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1014 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1015 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 1016 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1017 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1018 for (field = 0; field < numFields; ++field) { 1019 PetscObject obj; 1020 PetscClassId id; 1021 PetscInt Nc; 1022 1023 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 1024 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1025 if (id == PETSCFE_CLASSID) { 1026 PetscFE fe = (PetscFE) obj; 1027 1028 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1029 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1030 } else if (id == PETSCFV_CLASSID) { 1031 PetscFV fv = (PetscFV) obj; 1032 1033 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1034 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1035 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1036 numComponents += Nc; 1037 } 1038 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 1039 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 1040 ierr = PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);CHKERRQ(ierr); 1041 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1042 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1043 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1044 for (c = cStart; c < cEnd; ++c) { 1045 PetscScalar *x = NULL; 1046 PetscScalar elemDiff = 0.0; 1047 PetscInt qc = 0; 1048 1049 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);CHKERRQ(ierr); 1050 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1051 1052 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1053 PetscObject obj; 1054 PetscClassId id; 1055 void * const ctx = ctxs ? ctxs[field] : NULL; 1056 PetscInt Nb, Nc, q, fc; 1057 1058 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 1059 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1060 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1061 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1062 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1063 if (funcs[field]) { 1064 for (q = 0; q < Nq; ++q) { 1065 if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", (double)detJ[q], c, q); 1066 ierr = (*funcs[field])(coordDim, time, &coords[q*coordDim], Nc, funcVal, ctx); 1067 if (ierr) { 1068 PetscErrorCode ierr2; 1069 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 1070 ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2); 1071 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 1072 CHKERRQ(ierr); 1073 } 1074 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1075 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1076 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1077 for (fc = 0; fc < Nc; ++fc) { 1078 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 1079 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]; 1080 } 1081 } 1082 } 1083 fieldOffset += Nb; 1084 qc += Nc; 1085 } 1086 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1087 ierr = VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);CHKERRQ(ierr); 1088 } 1089 ierr = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr); 1090 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1091 ierr = VecSqrtAbs(D);CHKERRQ(ierr); 1092 PetscFunctionReturn(0); 1093 } 1094 1095 /*@C 1096 DMPlexComputeGradientClementInterpolant - This function computes the L2 projection of the cellwise gradient of a function u onto P1, and stores it in a Vec. 1097 1098 Input Parameters: 1099 + dm - The DM 1100 - LocX - The coefficient vector u_h 1101 1102 Output Parameter: 1103 . locC - A Vec which holds the Clement interpolant of the gradient 1104 1105 Notes: 1106 Add citation to (Clement, 1975) and definition of the interpolant 1107 \nabla u_h(v_i) = \sum_{T_i \in support(v_i)} |T_i| \nabla u_h(T_i) / \sum_{T_i \in support(v_i)} |T_i| where |T_i| is the cell volume 1108 1109 Level: developer 1110 1111 .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff() 1112 @*/ 1113 PetscErrorCode DMPlexComputeGradientClementInterpolant(DM dm, Vec locX, Vec locC) 1114 { 1115 DM_Plex *mesh = (DM_Plex *) dm->data; 1116 PetscInt debug = mesh->printFEM; 1117 DM dmC; 1118 PetscSection section; 1119 PetscQuadrature quad; 1120 PetscScalar *interpolant, *gradsum; 1121 PetscReal *coords, *detJ, *J, *invJ; 1122 const PetscReal *quadPoints, *quadWeights; 1123 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, vStart, vEnd, v, field, fieldOffset; 1124 PetscErrorCode ierr; 1125 1126 PetscFunctionBegin; 1127 ierr = VecGetDM(locC, &dmC);CHKERRQ(ierr); 1128 ierr = VecSet(locC, 0.0);CHKERRQ(ierr); 1129 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1130 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 1131 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 1132 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1133 for (field = 0; field < numFields; ++field) { 1134 PetscObject obj; 1135 PetscClassId id; 1136 PetscInt Nc; 1137 1138 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 1139 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1140 if (id == PETSCFE_CLASSID) { 1141 PetscFE fe = (PetscFE) obj; 1142 1143 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1144 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1145 } else if (id == PETSCFV_CLASSID) { 1146 PetscFV fv = (PetscFV) obj; 1147 1148 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1149 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1150 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1151 numComponents += Nc; 1152 } 1153 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 1154 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 1155 ierr = PetscMalloc6(coordDim*numComponents*2,&gradsum,coordDim*numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ);CHKERRQ(ierr); 1156 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1157 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1158 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1159 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1160 for (v = vStart; v < vEnd; ++v) { 1161 PetscScalar volsum = 0.0; 1162 PetscInt *star = NULL; 1163 PetscInt starSize, st, d, fc; 1164 1165 ierr = PetscMemzero(gradsum, coordDim*numComponents * sizeof(PetscScalar));CHKERRQ(ierr); 1166 ierr = DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1167 for (st = 0; st < starSize*2; st += 2) { 1168 const PetscInt cell = star[st]; 1169 PetscScalar *grad = &gradsum[coordDim*numComponents]; 1170 PetscScalar *x = NULL; 1171 PetscReal vol = 0.0; 1172 1173 if ((cell < cStart) || (cell >= cEnd)) continue; 1174 ierr = DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, J, invJ, detJ);CHKERRQ(ierr); 1175 ierr = DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr); 1176 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1177 PetscObject obj; 1178 PetscClassId id; 1179 PetscInt Nb, Nc, q, qc = 0; 1180 1181 ierr = PetscMemzero(grad, coordDim*numComponents * sizeof(PetscScalar));CHKERRQ(ierr); 1182 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 1183 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1184 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1185 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1186 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1187 for (q = 0; q < Nq; ++q) { 1188 if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", (double)detJ[q], cell, q); 1189 if (ierr) { 1190 PetscErrorCode ierr2; 1191 ierr2 = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr2); 1192 ierr2 = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr2); 1193 ierr2 = PetscFree4(interpolant,coords,detJ,J);CHKERRQ(ierr2); 1194 CHKERRQ(ierr); 1195 } 1196 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolateGradient_Static((PetscFE) obj, &x[fieldOffset], coordDim, invJ, NULL, q, interpolant);CHKERRQ(ierr);} 1197 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1198 for (fc = 0; fc < Nc; ++fc) { 1199 const PetscReal wt = quadWeights[q*qNc+qc+fc]; 1200 1201 for (d = 0; d < coordDim; ++d) grad[fc*coordDim+d] += interpolant[fc*dim+d]*wt*detJ[q]; 1202 } 1203 vol += quadWeights[q*qNc]*detJ[q]; 1204 } 1205 fieldOffset += Nb; 1206 qc += Nc; 1207 } 1208 ierr = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr); 1209 for (fc = 0; fc < numComponents; ++fc) { 1210 for (d = 0; d < coordDim; ++d) { 1211 gradsum[fc*coordDim+d] += grad[fc*coordDim+d]; 1212 } 1213 } 1214 volsum += vol; 1215 if (debug) { 1216 ierr = PetscPrintf(PETSC_COMM_SELF, "Cell %D gradient: [", cell);CHKERRQ(ierr); 1217 for (fc = 0; fc < numComponents; ++fc) { 1218 for (d = 0; d < coordDim; ++d) { 1219 if (fc || d > 0) {ierr = PetscPrintf(PETSC_COMM_SELF, ", ");CHKERRQ(ierr);} 1220 ierr = PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc*coordDim+d]));CHKERRQ(ierr); 1221 } 1222 } 1223 ierr = PetscPrintf(PETSC_COMM_SELF, "]\n");CHKERRQ(ierr); 1224 } 1225 } 1226 for (fc = 0; fc < numComponents; ++fc) { 1227 for (d = 0; d < coordDim; ++d) gradsum[fc*coordDim+d] /= volsum; 1228 } 1229 ierr = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1230 ierr = DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES);CHKERRQ(ierr); 1231 } 1232 ierr = PetscFree6(gradsum,interpolant,coords,detJ,J,invJ);CHKERRQ(ierr); 1233 PetscFunctionReturn(0); 1234 } 1235 1236 static PetscErrorCode DMPlexComputeIntegral_Internal(DM dm, Vec X, PetscInt cStart, PetscInt cEnd, PetscScalar *cintegral, void *user) 1237 { 1238 DM dmAux = NULL; 1239 PetscDS prob, probAux = NULL; 1240 PetscSection section, sectionAux; 1241 Vec locX, locA; 1242 PetscInt dim, numCells = cEnd - cStart, c, f; 1243 PetscBool useFVM = PETSC_FALSE; 1244 /* DS */ 1245 PetscInt Nf, totDim, *uOff, *uOff_x, numConstants; 1246 PetscInt NfAux, totDimAux, *aOff; 1247 PetscScalar *u, *a; 1248 const PetscScalar *constants; 1249 /* Geometry */ 1250 PetscFEGeom *cgeomFEM; 1251 DM dmGrad; 1252 PetscQuadrature affineQuad = NULL; 1253 Vec cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL; 1254 PetscFVCellGeom *cgeomFVM; 1255 const PetscScalar *lgrad; 1256 PetscInt maxDegree; 1257 DMField coordField; 1258 IS cellIS; 1259 PetscErrorCode ierr; 1260 1261 PetscFunctionBegin; 1262 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1263 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1264 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 1265 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1266 /* Determine which discretizations we have */ 1267 for (f = 0; f < Nf; ++f) { 1268 PetscObject obj; 1269 PetscClassId id; 1270 1271 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1272 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1273 if (id == PETSCFV_CLASSID) useFVM = PETSC_TRUE; 1274 } 1275 /* Get local solution with boundary values */ 1276 ierr = DMGetLocalVector(dm, &locX);CHKERRQ(ierr); 1277 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 1278 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);CHKERRQ(ierr); 1279 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);CHKERRQ(ierr); 1280 /* Read DS information */ 1281 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1282 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 1283 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 1284 ierr = ISCreateStride(PETSC_COMM_SELF,numCells,cStart,1,&cellIS);CHKERRQ(ierr); 1285 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 1286 /* Read Auxiliary DS information */ 1287 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1288 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);CHKERRQ(ierr); 1289 if (dmAux) { 1290 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1291 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 1292 ierr = DMGetSection(dmAux, §ionAux);CHKERRQ(ierr); 1293 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1294 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 1295 } 1296 /* Allocate data arrays */ 1297 ierr = PetscCalloc1(numCells*totDim, &u);CHKERRQ(ierr); 1298 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1299 /* Read out geometry */ 1300 ierr = DMGetCoordinateField(dm,&coordField);CHKERRQ(ierr); 1301 ierr = DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);CHKERRQ(ierr); 1302 if (maxDegree <= 1) { 1303 ierr = DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);CHKERRQ(ierr); 1304 if (affineQuad) { 1305 ierr = DMFieldCreateFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&cgeomFEM);CHKERRQ(ierr); 1306 } 1307 } 1308 if (useFVM) { 1309 PetscFV fv = NULL; 1310 Vec grad; 1311 PetscInt fStart, fEnd; 1312 PetscBool compGrad; 1313 1314 for (f = 0; f < Nf; ++f) { 1315 PetscObject obj; 1316 PetscClassId id; 1317 1318 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1319 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1320 if (id == PETSCFV_CLASSID) {fv = (PetscFV) obj; break;} 1321 } 1322 ierr = PetscFVGetComputeGradients(fv, &compGrad);CHKERRQ(ierr); 1323 ierr = PetscFVSetComputeGradients(fv, PETSC_TRUE);CHKERRQ(ierr); 1324 ierr = DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);CHKERRQ(ierr); 1325 ierr = DMPlexComputeGradientFVM(dm, fv, faceGeometryFVM, cellGeometryFVM, &dmGrad);CHKERRQ(ierr); 1326 ierr = PetscFVSetComputeGradients(fv, compGrad);CHKERRQ(ierr); 1327 ierr = VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1328 /* Reconstruct and limit cell gradients */ 1329 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1330 ierr = DMGetGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1331 ierr = DMPlexReconstructGradients_Internal(dm, fv, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);CHKERRQ(ierr); 1332 /* Communicate gradient values */ 1333 ierr = DMGetLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1334 ierr = DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1335 ierr = DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1336 ierr = DMRestoreGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1337 /* Handle non-essential (e.g. outflow) boundary values */ 1338 ierr = DMPlexInsertBoundaryValues(dm, PETSC_FALSE, locX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);CHKERRQ(ierr); 1339 ierr = VecGetArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1340 } 1341 /* Read out data from inputs */ 1342 for (c = cStart; c < cEnd; ++c) { 1343 PetscScalar *x = NULL; 1344 PetscInt i; 1345 1346 ierr = DMPlexVecGetClosure(dm, section, locX, c, NULL, &x);CHKERRQ(ierr); 1347 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 1348 ierr = DMPlexVecRestoreClosure(dm, section, locX, c, NULL, &x);CHKERRQ(ierr); 1349 if (dmAux) { 1350 ierr = DMPlexVecGetClosure(dmAux, sectionAux, locA, c, NULL, &x);CHKERRQ(ierr); 1351 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 1352 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, locA, c, NULL, &x);CHKERRQ(ierr); 1353 } 1354 } 1355 /* Do integration for each field */ 1356 for (f = 0; f < Nf; ++f) { 1357 PetscObject obj; 1358 PetscClassId id; 1359 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 1360 1361 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1362 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1363 if (id == PETSCFE_CLASSID) { 1364 PetscFE fe = (PetscFE) obj; 1365 PetscQuadrature q; 1366 PetscFEGeom *chunkGeom = NULL; 1367 PetscInt Nq, Nb; 1368 1369 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1370 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1371 ierr = PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 1372 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1373 blockSize = Nb*Nq; 1374 batchSize = numBlocks * blockSize; 1375 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1376 numChunks = numCells / (numBatches*batchSize); 1377 Ne = numChunks*numBatches*batchSize; 1378 Nr = numCells % (numBatches*batchSize); 1379 offset = numCells - Nr; 1380 if (!affineQuad) { 1381 ierr = DMFieldCreateFEGeom(coordField,cellIS,q,PETSC_FALSE,&cgeomFEM);CHKERRQ(ierr); 1382 } 1383 ierr = PetscFEGeomGetChunk(cgeomFEM,0,offset,&chunkGeom);CHKERRQ(ierr); 1384 ierr = PetscFEIntegrate(fe, prob, f, Ne, chunkGeom, u, probAux, a, cintegral);CHKERRQ(ierr); 1385 ierr = PetscFEGeomGetChunk(cgeomFEM,offset,numCells,&chunkGeom);CHKERRQ(ierr); 1386 ierr = PetscFEIntegrate(fe, prob, f, Nr, chunkGeom, &u[offset*totDim], probAux, &a[offset*totDimAux], &cintegral[offset*Nf]);CHKERRQ(ierr); 1387 ierr = PetscFEGeomRestoreChunk(cgeomFEM,offset,numCells,&chunkGeom);CHKERRQ(ierr); 1388 if (!affineQuad) { 1389 ierr = PetscFEGeomDestroy(&cgeomFEM);CHKERRQ(ierr); 1390 } 1391 } else if (id == PETSCFV_CLASSID) { 1392 PetscInt foff; 1393 PetscPointFunc obj_func; 1394 PetscScalar lint; 1395 1396 ierr = PetscDSGetObjective(prob, f, &obj_func);CHKERRQ(ierr); 1397 ierr = PetscDSGetFieldOffset(prob, f, &foff);CHKERRQ(ierr); 1398 if (obj_func) { 1399 for (c = 0; c < numCells; ++c) { 1400 PetscScalar *u_x; 1401 1402 ierr = DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);CHKERRQ(ierr); 1403 obj_func(dim, Nf, NfAux, uOff, uOff_x, &u[totDim*c+foff], NULL, u_x, aOff, NULL, &a[totDimAux*c], NULL, NULL, 0.0, cgeomFVM[c].centroid, numConstants, constants, &lint); 1404 cintegral[c*Nf+f] += PetscRealPart(lint)*cgeomFVM[c].volume; 1405 } 1406 } 1407 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 1408 } 1409 /* Cleanup data arrays */ 1410 if (useFVM) { 1411 ierr = VecRestoreArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1412 ierr = VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1413 ierr = DMRestoreLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1414 ierr = VecDestroy(&faceGeometryFVM);CHKERRQ(ierr); 1415 ierr = VecDestroy(&cellGeometryFVM);CHKERRQ(ierr); 1416 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr); 1417 } 1418 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1419 ierr = PetscFree(u);CHKERRQ(ierr); 1420 /* Cleanup */ 1421 if (affineQuad) { 1422 ierr = PetscFEGeomDestroy(&cgeomFEM);CHKERRQ(ierr); 1423 } 1424 ierr = PetscQuadratureDestroy(&affineQuad);CHKERRQ(ierr); 1425 ierr = ISDestroy(&cellIS);CHKERRQ(ierr); 1426 ierr = DMRestoreLocalVector(dm, &locX);CHKERRQ(ierr); 1427 PetscFunctionReturn(0); 1428 } 1429 1430 /*@ 1431 DMPlexComputeIntegralFEM - Form the integral over the domain from the global input X using pointwise functions specified by the user 1432 1433 Input Parameters: 1434 + dm - The mesh 1435 . X - Global input vector 1436 - user - The user context 1437 1438 Output Parameter: 1439 . integral - Integral for each field 1440 1441 Level: developer 1442 1443 .seealso: DMPlexComputeResidualFEM() 1444 @*/ 1445 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscScalar *integral, void *user) 1446 { 1447 DM_Plex *mesh = (DM_Plex *) dm->data; 1448 PetscScalar *cintegral, *lintegral; 1449 PetscInt Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell; 1450 PetscErrorCode ierr; 1451 1452 PetscFunctionBegin; 1453 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1454 PetscValidHeaderSpecific(X, VEC_CLASSID, 2); 1455 PetscValidPointer(integral, 3); 1456 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1457 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 1458 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr); 1459 ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr); 1460 ierr = DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);CHKERRQ(ierr); 1461 cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight]; 1462 /* TODO Introduce a loop over large chunks (right now this is a single chunk) */ 1463 ierr = PetscCalloc2(Nf, &lintegral, (cEnd-cStart)*Nf, &cintegral);CHKERRQ(ierr); 1464 ierr = DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);CHKERRQ(ierr); 1465 /* Sum up values */ 1466 for (cell = cStart; cell < cEnd; ++cell) { 1467 const PetscInt c = cell - cStart; 1468 1469 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);CHKERRQ(ierr);} 1470 for (f = 0; f < Nf; ++f) lintegral[f] += cintegral[c*Nf+f]; 1471 } 1472 ierr = MPIU_Allreduce(lintegral, integral, Nf, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr); 1473 if (mesh->printFEM) { 1474 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Integral:");CHKERRQ(ierr); 1475 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", (double) PetscRealPart(integral[f]));CHKERRQ(ierr);} 1476 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 1477 } 1478 ierr = PetscFree2(lintegral, cintegral);CHKERRQ(ierr); 1479 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1480 PetscFunctionReturn(0); 1481 } 1482 1483 /*@ 1484 DMPlexComputeCellwiseIntegralFEM - Form the vector of cellwise integrals F from the global input X using pointwise functions specified by the user 1485 1486 Input Parameters: 1487 + dm - The mesh 1488 . X - Global input vector 1489 - user - The user context 1490 1491 Output Parameter: 1492 . integral - Cellwise integrals for each field 1493 1494 Level: developer 1495 1496 .seealso: DMPlexComputeResidualFEM() 1497 @*/ 1498 PetscErrorCode DMPlexComputeCellwiseIntegralFEM(DM dm, Vec X, Vec F, void *user) 1499 { 1500 DM_Plex *mesh = (DM_Plex *) dm->data; 1501 DM dmF; 1502 PetscSection sectionF; 1503 PetscScalar *cintegral, *af; 1504 PetscInt Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell; 1505 PetscErrorCode ierr; 1506 1507 PetscFunctionBegin; 1508 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1509 PetscValidHeaderSpecific(X, VEC_CLASSID, 2); 1510 PetscValidHeaderSpecific(F, VEC_CLASSID, 3); 1511 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1512 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 1513 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr); 1514 ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr); 1515 ierr = DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);CHKERRQ(ierr); 1516 cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight]; 1517 /* TODO Introduce a loop over large chunks (right now this is a single chunk) */ 1518 ierr = PetscCalloc1((cEnd-cStart)*Nf, &cintegral);CHKERRQ(ierr); 1519 ierr = DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);CHKERRQ(ierr); 1520 /* Put values in F*/ 1521 ierr = VecGetDM(F, &dmF);CHKERRQ(ierr); 1522 ierr = DMGetSection(dmF, §ionF);CHKERRQ(ierr); 1523 ierr = VecGetArray(F, &af);CHKERRQ(ierr); 1524 for (cell = cStart; cell < cEnd; ++cell) { 1525 const PetscInt c = cell - cStart; 1526 PetscInt dof, off; 1527 1528 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);CHKERRQ(ierr);} 1529 ierr = PetscSectionGetDof(sectionF, cell, &dof);CHKERRQ(ierr); 1530 ierr = PetscSectionGetOffset(sectionF, cell, &off);CHKERRQ(ierr); 1531 if (dof != Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cell dofs %D != %D", dof, Nf); 1532 for (f = 0; f < Nf; ++f) af[off+f] = cintegral[c*Nf+f]; 1533 } 1534 ierr = VecRestoreArray(F, &af);CHKERRQ(ierr); 1535 ierr = PetscFree(cintegral);CHKERRQ(ierr); 1536 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1537 PetscFunctionReturn(0); 1538 } 1539 1540 static PetscErrorCode DMPlexComputeBdIntegral_Internal(DM dm, Vec locX, IS pointIS, 1541 void (*func)(PetscInt, PetscInt, PetscInt, 1542 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 1543 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 1544 PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 1545 PetscScalar *fintegral, void *user) 1546 { 1547 DM plex = NULL, plexA = NULL; 1548 PetscDS prob, probAux = NULL; 1549 PetscSection section, sectionAux = NULL; 1550 Vec locA = NULL; 1551 DMField coordField; 1552 PetscInt Nf, totDim, *uOff, *uOff_x; 1553 PetscInt NfAux = 0, totDimAux = 0, *aOff = NULL; 1554 PetscScalar *u, *a = NULL; 1555 const PetscScalar *constants; 1556 PetscInt numConstants, f; 1557 PetscErrorCode ierr; 1558 1559 PetscFunctionBegin; 1560 ierr = DMGetCoordinateField(dm, &coordField);CHKERRQ(ierr); 1561 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 1562 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1563 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1564 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1565 /* Determine which discretizations we have */ 1566 for (f = 0; f < Nf; ++f) { 1567 PetscObject obj; 1568 PetscClassId id; 1569 1570 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1571 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1572 if (id == PETSCFV_CLASSID) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Not supported for FVM (field %D)", f); 1573 } 1574 /* Read DS information */ 1575 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1576 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 1577 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 1578 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 1579 /* Read Auxiliary DS information */ 1580 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);CHKERRQ(ierr); 1581 if (locA) { 1582 DM dmAux; 1583 1584 ierr = VecGetDM(locA, &dmAux);CHKERRQ(ierr); 1585 ierr = DMConvert(dmAux, DMPLEX, &plexA);CHKERRQ(ierr); 1586 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1587 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 1588 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1589 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1590 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 1591 } 1592 /* Integrate over points */ 1593 { 1594 PetscFEGeom *fgeom, *chunkGeom = NULL; 1595 PetscInt maxDegree; 1596 PetscQuadrature qGeom = NULL; 1597 const PetscInt *points; 1598 PetscInt numFaces, face, Nq, field; 1599 PetscInt numChunks, chunkSize, chunk, Nr, offset; 1600 1601 ierr = ISGetLocalSize(pointIS, &numFaces);CHKERRQ(ierr); 1602 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 1603 ierr = PetscCalloc2(numFaces*totDim, &u, locA ? numFaces*totDimAux : 0, &a);CHKERRQ(ierr); 1604 ierr = DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree);CHKERRQ(ierr); 1605 for (field = 0; field < Nf; ++field) { 1606 PetscFE fe; 1607 1608 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 1609 if (maxDegree <= 1) {ierr = DMFieldCreateDefaultQuadrature(coordField, pointIS, &qGeom);CHKERRQ(ierr);} 1610 if (!qGeom) { 1611 ierr = PetscFEGetFaceQuadrature(fe, &qGeom);CHKERRQ(ierr); 1612 ierr = PetscObjectReference((PetscObject) qGeom);CHKERRQ(ierr); 1613 } 1614 ierr = PetscQuadratureGetData(qGeom, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 1615 ierr = DMPlexGetFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom);CHKERRQ(ierr); 1616 for (face = 0; face < numFaces; ++face) { 1617 const PetscInt point = points[face], *support, *cone; 1618 PetscScalar *x = NULL; 1619 PetscInt i, coneSize, faceLoc; 1620 1621 ierr = DMPlexGetSupport(dm, point, &support);CHKERRQ(ierr); 1622 ierr = DMPlexGetConeSize(dm, support[0], &coneSize);CHKERRQ(ierr); 1623 ierr = DMPlexGetCone(dm, support[0], &cone);CHKERRQ(ierr); 1624 for (faceLoc = 0; faceLoc < coneSize; ++faceLoc) if (cone[faceLoc] == point) break; 1625 if (faceLoc == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %D in cone of support[0] %D", face, support[0]); 1626 fgeom->face[face][0] = faceLoc; 1627 ierr = DMPlexVecGetClosure(plex, section, locX, support[0], NULL, &x);CHKERRQ(ierr); 1628 for (i = 0; i < totDim; ++i) u[face*totDim+i] = x[i]; 1629 ierr = DMPlexVecRestoreClosure(plex, section, locX, support[0], NULL, &x);CHKERRQ(ierr); 1630 if (locA) { 1631 PetscInt subp; 1632 ierr = DMPlexGetSubpoint(plexA, support[0], &subp);CHKERRQ(ierr); 1633 ierr = DMPlexVecGetClosure(plexA, sectionAux, locA, subp, NULL, &x);CHKERRQ(ierr); 1634 for (i = 0; i < totDimAux; ++i) a[f*totDimAux+i] = x[i]; 1635 ierr = DMPlexVecRestoreClosure(plexA, sectionAux, locA, subp, NULL, &x);CHKERRQ(ierr); 1636 } 1637 } 1638 /* Get blocking */ 1639 { 1640 PetscQuadrature q; 1641 PetscInt numBatches, batchSize, numBlocks, blockSize; 1642 PetscInt Nq, Nb; 1643 1644 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1645 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1646 ierr = PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 1647 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1648 blockSize = Nb*Nq; 1649 batchSize = numBlocks * blockSize; 1650 chunkSize = numBatches*batchSize; 1651 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1652 numChunks = numFaces / chunkSize; 1653 Nr = numFaces % chunkSize; 1654 offset = numFaces - Nr; 1655 } 1656 /* Do integration for each field */ 1657 for (chunk = 0; chunk < numChunks; ++chunk) { 1658 ierr = PetscFEGeomGetChunk(fgeom, chunk*chunkSize, (chunk+1)*chunkSize, &chunkGeom);CHKERRQ(ierr); 1659 ierr = PetscFEIntegrateBd(fe, prob, field, func, chunkSize, chunkGeom, u, probAux, a, fintegral);CHKERRQ(ierr); 1660 ierr = PetscFEGeomRestoreChunk(fgeom, 0, offset, &chunkGeom);CHKERRQ(ierr); 1661 } 1662 ierr = PetscFEGeomGetChunk(fgeom, offset, numFaces, &chunkGeom);CHKERRQ(ierr); 1663 ierr = PetscFEIntegrateBd(fe, prob, field, func, Nr, chunkGeom, &u[offset*totDim], probAux, a ? &a[offset*totDimAux] : NULL, &fintegral[offset*Nf]);CHKERRQ(ierr); 1664 ierr = PetscFEGeomRestoreChunk(fgeom, offset, numFaces, &chunkGeom);CHKERRQ(ierr); 1665 /* Cleanup data arrays */ 1666 ierr = DMPlexRestoreFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom);CHKERRQ(ierr); 1667 ierr = PetscQuadratureDestroy(&qGeom);CHKERRQ(ierr); 1668 ierr = PetscFree2(u, a);CHKERRQ(ierr); 1669 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 1670 } 1671 } 1672 if (plex) {ierr = DMDestroy(&plex);CHKERRQ(ierr);} 1673 if (plexA) {ierr = DMDestroy(&plexA);CHKERRQ(ierr);} 1674 PetscFunctionReturn(0); 1675 } 1676 1677 /*@ 1678 DMPlexComputeBdIntegral - Form the integral over the specified boundary from the global input X using pointwise functions specified by the user 1679 1680 Input Parameters: 1681 + dm - The mesh 1682 . X - Global input vector 1683 . label - The boundary DMLabel 1684 . numVals - The number of label values to use, or PETSC_DETERMINE for all values 1685 . vals - The label values to use, or PETSC_NULL for all values 1686 . func = The function to integrate along the boundary 1687 - user - The user context 1688 1689 Output Parameter: 1690 . integral - Integral for each field 1691 1692 Level: developer 1693 1694 .seealso: DMPlexComputeIntegralFEM(), DMPlexComputeBdResidualFEM() 1695 @*/ 1696 PetscErrorCode DMPlexComputeBdIntegral(DM dm, Vec X, DMLabel label, PetscInt numVals, const PetscInt vals[], 1697 void (*func)(PetscInt, PetscInt, PetscInt, 1698 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 1699 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 1700 PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 1701 PetscScalar *integral, void *user) 1702 { 1703 Vec locX; 1704 PetscSection section; 1705 DMLabel depthLabel; 1706 IS facetIS; 1707 PetscInt dim, Nf, f, v; 1708 PetscErrorCode ierr; 1709 1710 PetscFunctionBegin; 1711 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1712 PetscValidHeaderSpecific(X, VEC_CLASSID, 2); 1713 PetscValidPointer(label, 3); 1714 if (vals) PetscValidPointer(vals, 5); 1715 PetscValidPointer(integral, 6); 1716 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1717 ierr = DMPlexGetDepthLabel(dm, &depthLabel);CHKERRQ(ierr); 1718 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1719 ierr = DMLabelGetStratumIS(depthLabel, dim-1, &facetIS);CHKERRQ(ierr); 1720 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1721 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1722 /* Get local solution with boundary values */ 1723 ierr = DMGetLocalVector(dm, &locX);CHKERRQ(ierr); 1724 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 1725 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);CHKERRQ(ierr); 1726 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);CHKERRQ(ierr); 1727 /* Loop over label values */ 1728 ierr = PetscMemzero(integral, Nf * sizeof(PetscScalar));CHKERRQ(ierr); 1729 for (v = 0; v < numVals; ++v) { 1730 IS pointIS; 1731 PetscInt numFaces, face; 1732 PetscScalar *fintegral; 1733 1734 ierr = DMLabelGetStratumIS(label, vals[v], &pointIS);CHKERRQ(ierr); 1735 if (!pointIS) continue; /* No points with that id on this process */ 1736 { 1737 IS isectIS; 1738 1739 /* TODO: Special cases of ISIntersect where it is quick to check a priori if one is a superset of the other */ 1740 ierr = ISIntersect_Caching_Internal(facetIS, pointIS, &isectIS);CHKERRQ(ierr); 1741 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 1742 pointIS = isectIS; 1743 } 1744 ierr = ISGetLocalSize(pointIS, &numFaces);CHKERRQ(ierr); 1745 ierr = PetscCalloc1(numFaces*Nf, &fintegral);CHKERRQ(ierr); 1746 ierr = DMPlexComputeBdIntegral_Internal(dm, locX, pointIS, func, fintegral, user);CHKERRQ(ierr); 1747 /* Sum point contributions into integral */ 1748 for (f = 0; f < Nf; ++f) for (face = 0; face < numFaces; ++face) integral[f] += fintegral[face*Nf+f]; 1749 ierr = PetscFree(fintegral);CHKERRQ(ierr); 1750 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 1751 } 1752 ierr = DMRestoreLocalVector(dm, &locX);CHKERRQ(ierr); 1753 ierr = ISDestroy(&facetIS);CHKERRQ(ierr); 1754 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1755 PetscFunctionReturn(0); 1756 } 1757 1758 /*@ 1759 DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1760 1761 Input Parameters: 1762 + dmf - The fine mesh 1763 . dmc - The coarse mesh 1764 - user - The user context 1765 1766 Output Parameter: 1767 . In - The interpolation matrix 1768 1769 Level: developer 1770 1771 .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1772 @*/ 1773 PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user) 1774 { 1775 DM_Plex *mesh = (DM_Plex *) dmc->data; 1776 const char *name = "Interpolator"; 1777 PetscDS prob; 1778 PetscFE *feRef; 1779 PetscFV *fvRef; 1780 PetscSection fsection, fglobalSection; 1781 PetscSection csection, cglobalSection; 1782 PetscScalar *elemMat; 1783 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c; 1784 PetscInt cTotDim, rTotDim = 0; 1785 PetscErrorCode ierr; 1786 1787 PetscFunctionBegin; 1788 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1789 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1790 ierr = DMGetSection(dmf, &fsection);CHKERRQ(ierr); 1791 ierr = DMGetGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1792 ierr = DMGetSection(dmc, &csection);CHKERRQ(ierr); 1793 ierr = DMGetGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1794 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1795 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1796 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1797 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1798 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr); 1799 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1800 for (f = 0; f < Nf; ++f) { 1801 PetscObject obj; 1802 PetscClassId id; 1803 PetscInt rNb = 0, Nc = 0; 1804 1805 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1806 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1807 if (id == PETSCFE_CLASSID) { 1808 PetscFE fe = (PetscFE) obj; 1809 1810 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1811 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1812 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1813 } else if (id == PETSCFV_CLASSID) { 1814 PetscFV fv = (PetscFV) obj; 1815 PetscDualSpace Q; 1816 1817 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1818 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1819 ierr = PetscDualSpaceGetDimension(Q, &rNb);CHKERRQ(ierr); 1820 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1821 } 1822 rTotDim += rNb; 1823 } 1824 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1825 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 1826 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 1827 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1828 PetscDualSpace Qref; 1829 PetscQuadrature f; 1830 const PetscReal *qpoints, *qweights; 1831 PetscReal *points; 1832 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1833 1834 /* Compose points from all dual basis functionals */ 1835 if (feRef[fieldI]) { 1836 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1837 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 1838 } else { 1839 ierr = PetscFVGetDualSpace(fvRef[fieldI], &Qref);CHKERRQ(ierr); 1840 ierr = PetscFVGetNumComponents(fvRef[fieldI], &Nc);CHKERRQ(ierr); 1841 } 1842 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1843 for (i = 0; i < fpdim; ++i) { 1844 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1845 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1846 npoints += Np; 1847 } 1848 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1849 for (i = 0, k = 0; i < fpdim; ++i) { 1850 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1851 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1852 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1853 } 1854 1855 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1856 PetscObject obj; 1857 PetscClassId id; 1858 PetscReal *B; 1859 PetscInt NcJ = 0, cpdim = 0, j, qNc; 1860 1861 ierr = PetscDSGetDiscretization(prob, fieldJ, &obj);CHKERRQ(ierr); 1862 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1863 if (id == PETSCFE_CLASSID) { 1864 PetscFE fe = (PetscFE) obj; 1865 1866 /* Evaluate basis at points */ 1867 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 1868 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1869 /* For now, fields only interpolate themselves */ 1870 if (fieldI == fieldJ) { 1871 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); 1872 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1873 for (i = 0, k = 0; i < fpdim; ++i) { 1874 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1875 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);CHKERRQ(ierr); 1876 if (qNc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in quadrature %D does not match coarse field %D", qNc, NcJ); 1877 for (p = 0; p < Np; ++p, ++k) { 1878 for (j = 0; j < cpdim; ++j) { 1879 /* 1880 cTotDim: Total columns in element interpolation matrix, sum of number of dual basis functionals in each field 1881 offsetI, offsetJ: Offsets into the larger element interpolation matrix for different fields 1882 fpdim, i, cpdim, j: Dofs for fine and coarse grids, correspond to dual space basis functionals 1883 qNC, Nc, Ncj, c: Number of components in this field 1884 Np, p: Number of quad points in the fine grid functional i 1885 k: i*Np + p, overall point number for the interpolation 1886 */ 1887 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p*qNc+c]; 1888 } 1889 } 1890 } 1891 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1892 } 1893 } else if (id == PETSCFV_CLASSID) { 1894 PetscFV fv = (PetscFV) obj; 1895 1896 /* Evaluate constant function at points */ 1897 ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr); 1898 cpdim = 1; 1899 /* For now, fields only interpolate themselves */ 1900 if (fieldI == fieldJ) { 1901 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); 1902 for (i = 0, k = 0; i < fpdim; ++i) { 1903 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1904 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);CHKERRQ(ierr); 1905 if (qNc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in quadrature %D does not match coarse field %D", qNc, NcJ); 1906 for (p = 0; p < Np; ++p, ++k) { 1907 for (j = 0; j < cpdim; ++j) { 1908 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += 1.0*qweights[p*qNc+c]; 1909 } 1910 } 1911 } 1912 } 1913 } 1914 offsetJ += cpdim; 1915 } 1916 offsetI += fpdim; 1917 ierr = PetscFree(points);CHKERRQ(ierr); 1918 } 1919 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1920 /* Preallocate matrix */ 1921 { 1922 Mat preallocator; 1923 PetscScalar *vals; 1924 PetscInt *cellCIndices, *cellFIndices; 1925 PetscInt locRows, locCols, cell; 1926 1927 ierr = MatGetLocalSize(In, &locRows, &locCols);CHKERRQ(ierr); 1928 ierr = MatCreate(PetscObjectComm((PetscObject) In), &preallocator);CHKERRQ(ierr); 1929 ierr = MatSetType(preallocator, MATPREALLOCATOR);CHKERRQ(ierr); 1930 ierr = MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1931 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 1932 ierr = PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1933 for (cell = cStart; cell < cEnd; ++cell) { 1934 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1935 ierr = MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);CHKERRQ(ierr); 1936 } 1937 ierr = PetscFree3(vals,cellCIndices,cellFIndices);CHKERRQ(ierr); 1938 ierr = MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1939 ierr = MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1940 ierr = MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);CHKERRQ(ierr); 1941 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 1942 } 1943 /* Fill matrix */ 1944 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1945 for (c = cStart; c < cEnd; ++c) { 1946 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1947 } 1948 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1949 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1950 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1951 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1952 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1953 if (mesh->printFEM) { 1954 ierr = PetscPrintf(PetscObjectComm((PetscObject)In), "%s:\n", name);CHKERRQ(ierr); 1955 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1956 ierr = MatView(In, NULL);CHKERRQ(ierr); 1957 } 1958 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1959 PetscFunctionReturn(0); 1960 } 1961 1962 PetscErrorCode DMPlexComputeMassMatrixNested(DM dmc, DM dmf, Mat mass, void *user) 1963 { 1964 SETERRQ(PetscObjectComm((PetscObject) dmc), PETSC_ERR_SUP, "Laziness"); 1965 } 1966 1967 /*@ 1968 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 1969 1970 Input Parameters: 1971 + dmf - The fine mesh 1972 . dmc - The coarse mesh 1973 - user - The user context 1974 1975 Output Parameter: 1976 . In - The interpolation matrix 1977 1978 Level: developer 1979 1980 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1981 @*/ 1982 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 1983 { 1984 DM_Plex *mesh = (DM_Plex *) dmf->data; 1985 const char *name = "Interpolator"; 1986 PetscDS prob; 1987 PetscSection fsection, csection, globalFSection, globalCSection; 1988 PetscHSetIJ ht; 1989 PetscLayout rLayout; 1990 PetscInt *dnz, *onz; 1991 PetscInt locRows, rStart, rEnd; 1992 PetscReal *x, *v0, *J, *invJ, detJ; 1993 PetscReal *v0c, *Jc, *invJc, detJc; 1994 PetscScalar *elemMat; 1995 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1996 PetscErrorCode ierr; 1997 1998 PetscFunctionBegin; 1999 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2000 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 2001 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 2002 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 2003 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 2004 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 2005 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 2006 ierr = DMGetSection(dmf, &fsection);CHKERRQ(ierr); 2007 ierr = DMGetGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 2008 ierr = DMGetSection(dmc, &csection);CHKERRQ(ierr); 2009 ierr = DMGetGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 2010 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 2011 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 2012 ierr = PetscMalloc1(totDim, &elemMat);CHKERRQ(ierr); 2013 2014 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 2015 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 2016 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 2017 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 2018 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 2019 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 2020 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 2021 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 2022 ierr = PetscHSetIJCreate(&ht);CHKERRQ(ierr); 2023 for (field = 0; field < Nf; ++field) { 2024 PetscObject obj; 2025 PetscClassId id; 2026 PetscDualSpace Q = NULL; 2027 PetscQuadrature f; 2028 const PetscReal *qpoints; 2029 PetscInt Nc, Np, fpdim, i, d; 2030 2031 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 2032 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2033 if (id == PETSCFE_CLASSID) { 2034 PetscFE fe = (PetscFE) obj; 2035 2036 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 2037 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 2038 } else if (id == PETSCFV_CLASSID) { 2039 PetscFV fv = (PetscFV) obj; 2040 2041 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 2042 Nc = 1; 2043 } 2044 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 2045 /* For each fine grid cell */ 2046 for (cell = cStart; cell < cEnd; ++cell) { 2047 PetscInt *findices, *cindices; 2048 PetscInt numFIndices, numCIndices; 2049 2050 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2051 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 2052 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim); 2053 for (i = 0; i < fpdim; ++i) { 2054 Vec pointVec; 2055 PetscScalar *pV; 2056 PetscSF coarseCellSF = NULL; 2057 const PetscSFNode *coarseCells; 2058 PetscInt numCoarseCells, q, c; 2059 2060 /* Get points from the dual basis functional quadrature */ 2061 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 2062 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 2063 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 2064 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 2065 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2066 for (q = 0; q < Np; ++q) { 2067 const PetscReal xi0[3] = {-1., -1., -1.}; 2068 2069 /* Transform point to real space */ 2070 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x); 2071 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 2072 } 2073 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2074 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 2075 /* OPT: Pack all quad points from fine cell */ 2076 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 2077 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 2078 /* Update preallocation info */ 2079 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 2080 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 2081 { 2082 PetscHashIJKey key; 2083 PetscBool missing; 2084 2085 key.i = findices[i]; 2086 if (key.i >= 0) { 2087 /* Get indices for coarse elements */ 2088 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 2089 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2090 for (c = 0; c < numCIndices; ++c) { 2091 key.j = cindices[c]; 2092 if (key.j < 0) continue; 2093 ierr = PetscHSetIJQueryAdd(ht, key, &missing);CHKERRQ(ierr); 2094 if (missing) { 2095 if ((key.j >= rStart) && (key.j < rEnd)) ++dnz[key.i-rStart]; 2096 else ++onz[key.i-rStart]; 2097 } 2098 } 2099 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2100 } 2101 } 2102 } 2103 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 2104 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 2105 } 2106 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2107 } 2108 } 2109 ierr = PetscHSetIJDestroy(&ht);CHKERRQ(ierr); 2110 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 2111 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2112 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 2113 for (field = 0; field < Nf; ++field) { 2114 PetscObject obj; 2115 PetscClassId id; 2116 PetscDualSpace Q = NULL; 2117 PetscQuadrature f; 2118 const PetscReal *qpoints, *qweights; 2119 PetscInt Nc, qNc, Np, fpdim, i, d; 2120 2121 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 2122 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2123 if (id == PETSCFE_CLASSID) { 2124 PetscFE fe = (PetscFE) obj; 2125 2126 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 2127 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 2128 } else if (id == PETSCFV_CLASSID) { 2129 PetscFV fv = (PetscFV) obj; 2130 2131 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 2132 Nc = 1; 2133 } else SETERRQ1(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONG,"Unknown discretization type for field %d",field); 2134 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 2135 /* For each fine grid cell */ 2136 for (cell = cStart; cell < cEnd; ++cell) { 2137 PetscInt *findices, *cindices; 2138 PetscInt numFIndices, numCIndices; 2139 2140 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2141 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 2142 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim); 2143 for (i = 0; i < fpdim; ++i) { 2144 Vec pointVec; 2145 PetscScalar *pV; 2146 PetscSF coarseCellSF = NULL; 2147 const PetscSFNode *coarseCells; 2148 PetscInt numCoarseCells, cpdim, q, c, j; 2149 2150 /* Get points from the dual basis functional quadrature */ 2151 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 2152 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);CHKERRQ(ierr); 2153 if (qNc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in quadrature %D does not match coarse field %D", qNc, Nc); 2154 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 2155 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 2156 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2157 for (q = 0; q < Np; ++q) { 2158 const PetscReal xi0[3] = {-1., -1., -1.}; 2159 2160 /* Transform point to real space */ 2161 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x); 2162 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 2163 } 2164 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2165 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 2166 /* OPT: Read this out from preallocation information */ 2167 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 2168 /* Update preallocation info */ 2169 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 2170 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 2171 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2172 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 2173 PetscReal pVReal[3]; 2174 const PetscReal xi0[3] = {-1., -1., -1.}; 2175 2176 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2177 /* Transform points from real space to coarse reference space */ 2178 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 2179 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 2180 CoordinatesRealToRef(dim, dim, xi0, v0c, invJc, pVReal, x); 2181 2182 if (id == PETSCFE_CLASSID) { 2183 PetscFE fe = (PetscFE) obj; 2184 PetscReal *B; 2185 2186 /* Evaluate coarse basis on contained point */ 2187 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 2188 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 2189 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 2190 /* Get elemMat entries by multiplying by weight */ 2191 for (j = 0; j < cpdim; ++j) { 2192 for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*qweights[ccell*qNc + c]; 2193 } 2194 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 2195 } else { 2196 cpdim = 1; 2197 for (j = 0; j < cpdim; ++j) { 2198 for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c]; 2199 } 2200 } 2201 /* Update interpolator */ 2202 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 2203 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 2204 ierr = MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 2205 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2206 } 2207 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2208 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 2209 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 2210 } 2211 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2212 } 2213 } 2214 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 2215 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 2216 ierr = PetscFree(elemMat);CHKERRQ(ierr); 2217 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2218 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2219 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2220 PetscFunctionReturn(0); 2221 } 2222 2223 /*@ 2224 DMPlexComputeMassMatrixGeneral - Form the local portion of the mass matrix M from the coarse DM to a non-nested fine DM. 2225 2226 Input Parameters: 2227 + dmf - The fine mesh 2228 . dmc - The coarse mesh 2229 - user - The user context 2230 2231 Output Parameter: 2232 . mass - The mass matrix 2233 2234 Level: developer 2235 2236 .seealso: DMPlexComputeMassMatrixNested(), DMPlexComputeInterpolatorNested(), DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 2237 @*/ 2238 PetscErrorCode DMPlexComputeMassMatrixGeneral(DM dmc, DM dmf, Mat mass, void *user) 2239 { 2240 DM_Plex *mesh = (DM_Plex *) dmf->data; 2241 const char *name = "Mass Matrix"; 2242 PetscDS prob; 2243 PetscSection fsection, csection, globalFSection, globalCSection; 2244 PetscHSetIJ ht; 2245 PetscLayout rLayout; 2246 PetscInt *dnz, *onz; 2247 PetscInt locRows, rStart, rEnd; 2248 PetscReal *x, *v0, *J, *invJ, detJ; 2249 PetscReal *v0c, *Jc, *invJc, detJc; 2250 PetscScalar *elemMat; 2251 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 2252 PetscErrorCode ierr; 2253 2254 PetscFunctionBegin; 2255 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 2256 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 2257 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 2258 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 2259 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 2260 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 2261 ierr = DMGetSection(dmf, &fsection);CHKERRQ(ierr); 2262 ierr = DMGetGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 2263 ierr = DMGetSection(dmc, &csection);CHKERRQ(ierr); 2264 ierr = DMGetGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 2265 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 2266 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 2267 ierr = PetscMalloc1(totDim, &elemMat);CHKERRQ(ierr); 2268 2269 ierr = MatGetLocalSize(mass, &locRows, NULL);CHKERRQ(ierr); 2270 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) mass), &rLayout);CHKERRQ(ierr); 2271 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 2272 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 2273 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 2274 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 2275 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 2276 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 2277 ierr = PetscHSetIJCreate(&ht);CHKERRQ(ierr); 2278 for (field = 0; field < Nf; ++field) { 2279 PetscObject obj; 2280 PetscClassId id; 2281 PetscQuadrature quad; 2282 const PetscReal *qpoints; 2283 PetscInt Nq, Nc, i, d; 2284 2285 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 2286 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2287 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetQuadrature((PetscFE) obj, &quad);CHKERRQ(ierr);} 2288 else {ierr = PetscFVGetQuadrature((PetscFV) obj, &quad);CHKERRQ(ierr);} 2289 ierr = PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, NULL);CHKERRQ(ierr); 2290 /* For each fine grid cell */ 2291 for (cell = cStart; cell < cEnd; ++cell) { 2292 Vec pointVec; 2293 PetscScalar *pV; 2294 PetscSF coarseCellSF = NULL; 2295 const PetscSFNode *coarseCells; 2296 PetscInt numCoarseCells, q, c; 2297 PetscInt *findices, *cindices; 2298 PetscInt numFIndices, numCIndices; 2299 2300 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2301 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 2302 /* Get points from the quadrature */ 2303 ierr = VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);CHKERRQ(ierr); 2304 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 2305 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2306 for (q = 0; q < Nq; ++q) { 2307 const PetscReal xi0[3] = {-1., -1., -1.}; 2308 2309 /* Transform point to real space */ 2310 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x); 2311 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 2312 } 2313 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2314 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 2315 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 2316 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 2317 /* Update preallocation info */ 2318 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 2319 if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 2320 { 2321 PetscHashIJKey key; 2322 PetscBool missing; 2323 2324 for (i = 0; i < numFIndices; ++i) { 2325 key.i = findices[i]; 2326 if (key.i >= 0) { 2327 /* Get indices for coarse elements */ 2328 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 2329 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2330 for (c = 0; c < numCIndices; ++c) { 2331 key.j = cindices[c]; 2332 if (key.j < 0) continue; 2333 ierr = PetscHSetIJQueryAdd(ht, key, &missing);CHKERRQ(ierr); 2334 if (missing) { 2335 if ((key.j >= rStart) && (key.j < rEnd)) ++dnz[key.i-rStart]; 2336 else ++onz[key.i-rStart]; 2337 } 2338 } 2339 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2340 } 2341 } 2342 } 2343 } 2344 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 2345 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 2346 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2347 } 2348 } 2349 ierr = PetscHSetIJDestroy(&ht);CHKERRQ(ierr); 2350 ierr = MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 2351 ierr = MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2352 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 2353 for (field = 0; field < Nf; ++field) { 2354 PetscObject obj; 2355 PetscClassId id; 2356 PetscQuadrature quad; 2357 PetscReal *Bfine; 2358 const PetscReal *qpoints, *qweights; 2359 PetscInt Nq, Nc, i, d; 2360 2361 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 2362 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2363 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetQuadrature((PetscFE) obj, &quad);CHKERRQ(ierr);ierr = PetscFEGetDefaultTabulation((PetscFE) obj, &Bfine, NULL, NULL);CHKERRQ(ierr);} 2364 else {ierr = PetscFVGetQuadrature((PetscFV) obj, &quad);CHKERRQ(ierr);} 2365 ierr = PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, &qweights);CHKERRQ(ierr); 2366 /* For each fine grid cell */ 2367 for (cell = cStart; cell < cEnd; ++cell) { 2368 Vec pointVec; 2369 PetscScalar *pV; 2370 PetscSF coarseCellSF = NULL; 2371 const PetscSFNode *coarseCells; 2372 PetscInt numCoarseCells, cpdim, q, c, j; 2373 PetscInt *findices, *cindices; 2374 PetscInt numFIndices, numCIndices; 2375 2376 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2377 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 2378 /* Get points from the quadrature */ 2379 ierr = VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);CHKERRQ(ierr); 2380 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 2381 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2382 for (q = 0; q < Nq; ++q) { 2383 const PetscReal xi0[3] = {-1., -1., -1.}; 2384 2385 /* Transform point to real space */ 2386 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x); 2387 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 2388 } 2389 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2390 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 2391 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 2392 /* Update matrix */ 2393 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 2394 if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 2395 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2396 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 2397 PetscReal pVReal[3]; 2398 const PetscReal xi0[3] = {-1., -1., -1.}; 2399 2400 2401 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2402 /* Transform points from real space to coarse reference space */ 2403 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 2404 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 2405 CoordinatesRealToRef(dim, dim, xi0, v0c, invJc, pVReal, x); 2406 2407 if (id == PETSCFE_CLASSID) { 2408 PetscFE fe = (PetscFE) obj; 2409 PetscReal *B; 2410 2411 /* Evaluate coarse basis on contained point */ 2412 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 2413 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 2414 /* Get elemMat entries by multiplying by weight */ 2415 for (i = 0; i < numFIndices; ++i) { 2416 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 2417 for (j = 0; j < cpdim; ++j) { 2418 for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*Bfine[(ccell*numFIndices + i)*Nc + c]*qweights[ccell*Nc + c]*detJ; 2419 } 2420 /* Update interpolator */ 2421 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 2422 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 2423 ierr = MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);CHKERRQ(ierr); 2424 } 2425 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 2426 } else { 2427 cpdim = 1; 2428 for (i = 0; i < numFIndices; ++i) { 2429 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 2430 for (j = 0; j < cpdim; ++j) { 2431 for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*1.0*qweights[ccell*Nc + c]*detJ; 2432 } 2433 /* Update interpolator */ 2434 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 2435 ierr = PetscPrintf(PETSC_COMM_SELF, "Nq: %d %d Nf: %d %d Nc: %d %d\n", ccell, Nq, i, numFIndices, j, numCIndices);CHKERRQ(ierr); 2436 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 2437 ierr = MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);CHKERRQ(ierr); 2438 } 2439 } 2440 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2441 } 2442 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2443 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 2444 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 2445 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2446 } 2447 } 2448 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 2449 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 2450 ierr = PetscFree(elemMat);CHKERRQ(ierr); 2451 ierr = MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2452 ierr = MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2453 PetscFunctionReturn(0); 2454 } 2455 2456 /*@ 2457 DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns 2458 2459 Input Parameters: 2460 + dmc - The coarse mesh 2461 - dmf - The fine mesh 2462 - user - The user context 2463 2464 Output Parameter: 2465 . sc - The mapping 2466 2467 Level: developer 2468 2469 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 2470 @*/ 2471 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 2472 { 2473 PetscDS prob; 2474 PetscFE *feRef; 2475 PetscFV *fvRef; 2476 Vec fv, cv; 2477 IS fis, cis; 2478 PetscSection fsection, fglobalSection, csection, cglobalSection; 2479 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 2480 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m; 2481 PetscBool *needAvg; 2482 PetscErrorCode ierr; 2483 2484 PetscFunctionBegin; 2485 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2486 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 2487 ierr = DMGetSection(dmf, &fsection);CHKERRQ(ierr); 2488 ierr = DMGetGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 2489 ierr = DMGetSection(dmc, &csection);CHKERRQ(ierr); 2490 ierr = DMGetGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 2491 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 2492 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 2493 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 2494 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 2495 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 2496 ierr = PetscCalloc3(Nf,&feRef,Nf,&fvRef,Nf,&needAvg);CHKERRQ(ierr); 2497 for (f = 0; f < Nf; ++f) { 2498 PetscObject obj; 2499 PetscClassId id; 2500 PetscInt fNb = 0, Nc = 0; 2501 2502 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 2503 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2504 if (id == PETSCFE_CLASSID) { 2505 PetscFE fe = (PetscFE) obj; 2506 PetscSpace sp; 2507 PetscInt maxDegree; 2508 2509 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 2510 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 2511 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 2512 ierr = PetscFEGetBasisSpace(fe, &sp);CHKERRQ(ierr); 2513 ierr = PetscSpaceGetDegree(sp, NULL, &maxDegree);CHKERRQ(ierr); 2514 if (!maxDegree) needAvg[f] = PETSC_TRUE; 2515 } else if (id == PETSCFV_CLASSID) { 2516 PetscFV fv = (PetscFV) obj; 2517 PetscDualSpace Q; 2518 2519 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 2520 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 2521 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 2522 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 2523 needAvg[f] = PETSC_TRUE; 2524 } 2525 fTotDim += fNb; 2526 } 2527 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 2528 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 2529 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 2530 PetscFE feC; 2531 PetscFV fvC; 2532 PetscDualSpace QF, QC; 2533 PetscInt order = -1, NcF, NcC, fpdim, cpdim; 2534 2535 if (feRef[field]) { 2536 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 2537 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 2538 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 2539 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 2540 ierr = PetscDualSpaceGetOrder(QF, &order);CHKERRQ(ierr); 2541 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 2542 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 2543 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 2544 } else { 2545 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 2546 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 2547 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 2548 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 2549 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 2550 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 2551 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 2552 } 2553 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); 2554 for (c = 0; c < cpdim; ++c) { 2555 PetscQuadrature cfunc; 2556 const PetscReal *cqpoints, *cqweights; 2557 PetscInt NqcC, NpC; 2558 PetscBool found = PETSC_FALSE; 2559 2560 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 2561 ierr = PetscQuadratureGetData(cfunc, NULL, &NqcC, &NpC, &cqpoints, &cqweights);CHKERRQ(ierr); 2562 if (NqcC != NcC) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of quadrature components %D must match number of field components", NqcC, NcC); 2563 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 2564 for (f = 0; f < fpdim; ++f) { 2565 PetscQuadrature ffunc; 2566 const PetscReal *fqpoints, *fqweights; 2567 PetscReal sum = 0.0; 2568 PetscInt NqcF, NpF; 2569 2570 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 2571 ierr = PetscQuadratureGetData(ffunc, NULL, &NqcF, &NpF, &fqpoints, &fqweights);CHKERRQ(ierr); 2572 if (NqcF != NcF) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of quadrature components %D must match number of field components", NqcF, NcF); 2573 if (NpC != NpF) continue; 2574 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 2575 if (sum > 1.0e-9) continue; 2576 for (d = 0; d < NcC; ++d) sum += PetscAbsReal(cqweights[d]*fqweights[d]); 2577 if (sum < 1.0e-9) continue; 2578 cmap[offsetC+c] = offsetF+f; 2579 found = PETSC_TRUE; 2580 break; 2581 } 2582 if (!found) { 2583 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 2584 if (fvRef[field] || (feRef[field] && order == 0)) { 2585 cmap[offsetC+c] = offsetF+0; 2586 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 2587 } 2588 } 2589 offsetC += cpdim; 2590 offsetF += fpdim; 2591 } 2592 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 2593 ierr = PetscFree3(feRef,fvRef,needAvg);CHKERRQ(ierr); 2594 2595 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 2596 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 2597 ierr = VecGetOwnershipRange(cv, &startC, &endC);CHKERRQ(ierr); 2598 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 2599 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 2600 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 2601 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 2602 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 2603 for (c = cStart; c < cEnd; ++c) { 2604 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 2605 for (d = 0; d < cTotDim; ++d) { 2606 if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue; 2607 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]]); 2608 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 2609 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 2610 } 2611 } 2612 ierr = PetscFree(cmap);CHKERRQ(ierr); 2613 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 2614 2615 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 2616 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 2617 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 2618 ierr = ISDestroy(&cis);CHKERRQ(ierr); 2619 ierr = ISDestroy(&fis);CHKERRQ(ierr); 2620 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 2621 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 2622 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2623 PetscFunctionReturn(0); 2624 } 2625 2626 /*@C 2627 DMPlexGetCellFields - Retrieve the field values values for a chunk of cells 2628 2629 Input Parameters: 2630 + dm - The DM 2631 . cellIS - The cells to include 2632 . locX - A local vector with the solution fields 2633 . locX_t - A local vector with solution field time derivatives, or NULL 2634 - locA - A local vector with auxiliary fields, or NULL 2635 2636 Output Parameters: 2637 + u - The field coefficients 2638 . u_t - The fields derivative coefficients 2639 - a - The auxiliary field coefficients 2640 2641 Level: developer 2642 2643 .seealso: DMPlexGetFaceFields() 2644 @*/ 2645 PetscErrorCode DMPlexGetCellFields(DM dm, IS cellIS, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a) 2646 { 2647 DM plex, plexA = NULL; 2648 PetscSection section, sectionAux; 2649 PetscDS prob; 2650 const PetscInt *cells; 2651 PetscInt cStart, cEnd, numCells, totDim, totDimAux, c; 2652 PetscErrorCode ierr; 2653 2654 PetscFunctionBegin; 2655 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2656 PetscValidHeaderSpecific(locX, VEC_CLASSID, 4); 2657 if (locX_t) {PetscValidHeaderSpecific(locX_t, VEC_CLASSID, 5);} 2658 if (locA) {PetscValidHeaderSpecific(locA, VEC_CLASSID, 6);} 2659 PetscValidPointer(u, 7); 2660 PetscValidPointer(u_t, 8); 2661 PetscValidPointer(a, 9); 2662 ierr = DMPlexConvertPlex(dm, &plex, PETSC_FALSE);CHKERRQ(ierr); 2663 ierr = ISGetPointRange(cellIS, &cStart, &cEnd, &cells);CHKERRQ(ierr); 2664 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 2665 ierr = DMGetCellDS(dm, cStart, &prob);CHKERRQ(ierr); 2666 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 2667 if (locA) { 2668 DM dmAux; 2669 PetscDS probAux; 2670 2671 ierr = VecGetDM(locA, &dmAux);CHKERRQ(ierr); 2672 ierr = DMPlexConvertPlex(dmAux, &plexA, PETSC_FALSE);CHKERRQ(ierr); 2673 ierr = DMGetSection(dmAux, §ionAux);CHKERRQ(ierr); 2674 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 2675 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 2676 } 2677 numCells = cEnd - cStart; 2678 ierr = DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, u);CHKERRQ(ierr); 2679 if (locX_t) {ierr = DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, u_t);CHKERRQ(ierr);} else {*u_t = NULL;} 2680 if (locA) {ierr = DMGetWorkArray(dm, numCells*totDimAux, MPIU_SCALAR, a);CHKERRQ(ierr);} else {*a = NULL;} 2681 for (c = cStart; c < cEnd; ++c) { 2682 const PetscInt cell = cells ? cells[c] : c; 2683 const PetscInt cind = c - cStart; 2684 PetscScalar *x = NULL, *x_t = NULL, *ul = *u, *ul_t = *u_t, *al = *a; 2685 PetscInt i; 2686 2687 ierr = DMPlexVecGetClosure(plex, section, locX, cell, NULL, &x);CHKERRQ(ierr); 2688 for (i = 0; i < totDim; ++i) ul[cind*totDim+i] = x[i]; 2689 ierr = DMPlexVecRestoreClosure(plex, section, locX, cell, NULL, &x);CHKERRQ(ierr); 2690 if (locX_t) { 2691 ierr = DMPlexVecGetClosure(plex, section, locX_t, cell, NULL, &x_t);CHKERRQ(ierr); 2692 for (i = 0; i < totDim; ++i) ul_t[cind*totDim+i] = x_t[i]; 2693 ierr = DMPlexVecRestoreClosure(plex, section, locX_t, cell, NULL, &x_t);CHKERRQ(ierr); 2694 } 2695 if (locA) { 2696 PetscInt subcell; 2697 ierr = DMPlexGetAuxiliaryPoint(plex, plexA, cell, &subcell);CHKERRQ(ierr); 2698 ierr = DMPlexVecGetClosure(plexA, sectionAux, locA, subcell, NULL, &x);CHKERRQ(ierr); 2699 for (i = 0; i < totDimAux; ++i) al[cind*totDimAux+i] = x[i]; 2700 ierr = DMPlexVecRestoreClosure(plexA, sectionAux, locA, subcell, NULL, &x);CHKERRQ(ierr); 2701 } 2702 } 2703 ierr = DMDestroy(&plex);CHKERRQ(ierr); 2704 if (locA) {ierr = DMDestroy(&plexA);CHKERRQ(ierr);} 2705 ierr = ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);CHKERRQ(ierr); 2706 PetscFunctionReturn(0); 2707 } 2708 2709 /*@C 2710 DMPlexRestoreCellFields - Restore the field values values for a chunk of cells 2711 2712 Input Parameters: 2713 + dm - The DM 2714 . cellIS - The cells to include 2715 . locX - A local vector with the solution fields 2716 . locX_t - A local vector with solution field time derivatives, or NULL 2717 - locA - A local vector with auxiliary fields, or NULL 2718 2719 Output Parameters: 2720 + u - The field coefficients 2721 . u_t - The fields derivative coefficients 2722 - a - The auxiliary field coefficients 2723 2724 Level: developer 2725 2726 .seealso: DMPlexGetFaceFields() 2727 @*/ 2728 PetscErrorCode DMPlexRestoreCellFields(DM dm, IS cellIS, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a) 2729 { 2730 PetscErrorCode ierr; 2731 2732 PetscFunctionBegin; 2733 ierr = DMRestoreWorkArray(dm, 0, MPIU_SCALAR, u);CHKERRQ(ierr); 2734 if (locX_t) {ierr = DMRestoreWorkArray(dm, 0, MPIU_SCALAR, u_t);CHKERRQ(ierr);} 2735 if (locA) {ierr = DMRestoreWorkArray(dm, 0, MPIU_SCALAR, a);CHKERRQ(ierr);} 2736 PetscFunctionReturn(0); 2737 } 2738 2739 /*@C 2740 DMPlexGetFaceFields - Retrieve the field values values for a chunk of faces 2741 2742 Input Parameters: 2743 + dm - The DM 2744 . fStart - The first face to include 2745 . fEnd - The first face to exclude 2746 . locX - A local vector with the solution fields 2747 . locX_t - A local vector with solution field time derivatives, or NULL 2748 . faceGeometry - A local vector with face geometry 2749 . cellGeometry - A local vector with cell geometry 2750 - locaGrad - A local vector with field gradients, or NULL 2751 2752 Output Parameters: 2753 + Nface - The number of faces with field values 2754 . uL - The field values at the left side of the face 2755 - uR - The field values at the right side of the face 2756 2757 Level: developer 2758 2759 .seealso: DMPlexGetCellFields() 2760 @*/ 2761 PetscErrorCode DMPlexGetFaceFields(DM dm, PetscInt fStart, PetscInt fEnd, Vec locX, Vec locX_t, Vec faceGeometry, Vec cellGeometry, Vec locGrad, PetscInt *Nface, PetscScalar **uL, PetscScalar **uR) 2762 { 2763 DM dmFace, dmCell, dmGrad = NULL; 2764 PetscSection section; 2765 PetscDS prob; 2766 DMLabel ghostLabel; 2767 const PetscScalar *facegeom, *cellgeom, *x, *lgrad; 2768 PetscBool *isFE; 2769 PetscInt dim, Nf, f, Nc, numFaces = fEnd - fStart, iface, face; 2770 PetscErrorCode ierr; 2771 2772 PetscFunctionBegin; 2773 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2774 PetscValidHeaderSpecific(locX, VEC_CLASSID, 4); 2775 if (locX_t) {PetscValidHeaderSpecific(locX_t, VEC_CLASSID, 5);} 2776 PetscValidHeaderSpecific(faceGeometry, VEC_CLASSID, 6); 2777 PetscValidHeaderSpecific(cellGeometry, VEC_CLASSID, 7); 2778 if (locGrad) {PetscValidHeaderSpecific(locGrad, VEC_CLASSID, 8);} 2779 PetscValidPointer(uL, 9); 2780 PetscValidPointer(uR, 10); 2781 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2782 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 2783 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 2784 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 2785 ierr = PetscDSGetTotalComponents(prob, &Nc);CHKERRQ(ierr); 2786 ierr = PetscMalloc1(Nf, &isFE);CHKERRQ(ierr); 2787 for (f = 0; f < Nf; ++f) { 2788 PetscObject obj; 2789 PetscClassId id; 2790 2791 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 2792 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2793 if (id == PETSCFE_CLASSID) {isFE[f] = PETSC_TRUE;} 2794 else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;} 2795 else {isFE[f] = PETSC_FALSE;} 2796 } 2797 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 2798 ierr = VecGetArrayRead(locX, &x);CHKERRQ(ierr); 2799 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 2800 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 2801 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 2802 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 2803 if (locGrad) { 2804 ierr = VecGetDM(locGrad, &dmGrad);CHKERRQ(ierr); 2805 ierr = VecGetArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 2806 } 2807 ierr = DMGetWorkArray(dm, numFaces*Nc, MPIU_SCALAR, uL);CHKERRQ(ierr); 2808 ierr = DMGetWorkArray(dm, numFaces*Nc, MPIU_SCALAR, uR);CHKERRQ(ierr); 2809 /* Right now just eat the extra work for FE (could make a cell loop) */ 2810 for (face = fStart, iface = 0; face < fEnd; ++face) { 2811 const PetscInt *cells; 2812 PetscFVFaceGeom *fg; 2813 PetscFVCellGeom *cgL, *cgR; 2814 PetscScalar *xL, *xR, *gL, *gR; 2815 PetscScalar *uLl = *uL, *uRl = *uR; 2816 PetscInt ghost, nsupp, nchild; 2817 2818 ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr); 2819 ierr = DMPlexGetSupportSize(dm, face, &nsupp);CHKERRQ(ierr); 2820 ierr = DMPlexGetTreeChildren(dm, face, &nchild, NULL);CHKERRQ(ierr); 2821 if (ghost >= 0 || nsupp > 2 || nchild > 0) continue; 2822 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 2823 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 2824 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);CHKERRQ(ierr); 2825 ierr = DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);CHKERRQ(ierr); 2826 for (f = 0; f < Nf; ++f) { 2827 PetscInt off; 2828 2829 ierr = PetscDSGetComponentOffset(prob, f, &off);CHKERRQ(ierr); 2830 if (isFE[f]) { 2831 const PetscInt *cone; 2832 PetscInt comp, coneSizeL, coneSizeR, faceLocL, faceLocR, ldof, rdof, d; 2833 2834 xL = xR = NULL; 2835 ierr = PetscSectionGetFieldComponents(section, f, &comp);CHKERRQ(ierr); 2836 ierr = DMPlexVecGetClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);CHKERRQ(ierr); 2837 ierr = DMPlexVecGetClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);CHKERRQ(ierr); 2838 ierr = DMPlexGetCone(dm, cells[0], &cone);CHKERRQ(ierr); 2839 ierr = DMPlexGetConeSize(dm, cells[0], &coneSizeL);CHKERRQ(ierr); 2840 for (faceLocL = 0; faceLocL < coneSizeL; ++faceLocL) if (cone[faceLocL] == face) break; 2841 ierr = DMPlexGetCone(dm, cells[1], &cone);CHKERRQ(ierr); 2842 ierr = DMPlexGetConeSize(dm, cells[1], &coneSizeR);CHKERRQ(ierr); 2843 for (faceLocR = 0; faceLocR < coneSizeR; ++faceLocR) if (cone[faceLocR] == face) break; 2844 if (faceLocL == coneSizeL && faceLocR == coneSizeR) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %d in cone of cell %d or cell %d", face, cells[0], cells[1]); 2845 /* Check that FEM field has values in the right cell (sometimes its an FV ghost cell) */ 2846 /* TODO: this is a hack that might not be right for nonconforming */ 2847 if (faceLocL < coneSizeL) { 2848 ierr = EvaluateFaceFields(prob, f, faceLocL, xL, &uLl[iface*Nc+off]);CHKERRQ(ierr); 2849 if (rdof == ldof && faceLocR < coneSizeR) {ierr = EvaluateFaceFields(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);CHKERRQ(ierr);} 2850 else {for(d = 0; d < comp; ++d) uRl[iface*Nc+off+d] = uLl[iface*Nc+off+d];} 2851 } 2852 else { 2853 ierr = EvaluateFaceFields(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);CHKERRQ(ierr); 2854 ierr = PetscSectionGetFieldComponents(section, f, &comp);CHKERRQ(ierr); 2855 for(d = 0; d < comp; ++d) uLl[iface*Nc+off+d] = uRl[iface*Nc+off+d]; 2856 } 2857 ierr = DMPlexVecRestoreClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);CHKERRQ(ierr); 2858 ierr = DMPlexVecRestoreClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);CHKERRQ(ierr); 2859 } else { 2860 PetscFV fv; 2861 PetscInt numComp, c; 2862 2863 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fv);CHKERRQ(ierr); 2864 ierr = PetscFVGetNumComponents(fv, &numComp);CHKERRQ(ierr); 2865 ierr = DMPlexPointLocalFieldRead(dm, cells[0], f, x, &xL);CHKERRQ(ierr); 2866 ierr = DMPlexPointLocalFieldRead(dm, cells[1], f, x, &xR);CHKERRQ(ierr); 2867 if (dmGrad) { 2868 PetscReal dxL[3], dxR[3]; 2869 2870 ierr = DMPlexPointLocalRead(dmGrad, cells[0], lgrad, &gL);CHKERRQ(ierr); 2871 ierr = DMPlexPointLocalRead(dmGrad, cells[1], lgrad, &gR);CHKERRQ(ierr); 2872 DMPlex_WaxpyD_Internal(dim, -1, cgL->centroid, fg->centroid, dxL); 2873 DMPlex_WaxpyD_Internal(dim, -1, cgR->centroid, fg->centroid, dxR); 2874 for (c = 0; c < numComp; ++c) { 2875 uLl[iface*Nc+off+c] = xL[c] + DMPlex_DotD_Internal(dim, &gL[c*dim], dxL); 2876 uRl[iface*Nc+off+c] = xR[c] + DMPlex_DotD_Internal(dim, &gR[c*dim], dxR); 2877 } 2878 } else { 2879 for (c = 0; c < numComp; ++c) { 2880 uLl[iface*Nc+off+c] = xL[c]; 2881 uRl[iface*Nc+off+c] = xR[c]; 2882 } 2883 } 2884 } 2885 } 2886 ++iface; 2887 } 2888 *Nface = iface; 2889 ierr = VecRestoreArrayRead(locX, &x);CHKERRQ(ierr); 2890 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 2891 ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 2892 if (locGrad) { 2893 ierr = VecRestoreArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 2894 } 2895 ierr = PetscFree(isFE);CHKERRQ(ierr); 2896 PetscFunctionReturn(0); 2897 } 2898 2899 /*@C 2900 DMPlexRestoreFaceFields - Restore the field values values for a chunk of faces 2901 2902 Input Parameters: 2903 + dm - The DM 2904 . fStart - The first face to include 2905 . fEnd - The first face to exclude 2906 . locX - A local vector with the solution fields 2907 . locX_t - A local vector with solution field time derivatives, or NULL 2908 . faceGeometry - A local vector with face geometry 2909 . cellGeometry - A local vector with cell geometry 2910 - locaGrad - A local vector with field gradients, or NULL 2911 2912 Output Parameters: 2913 + Nface - The number of faces with field values 2914 . uL - The field values at the left side of the face 2915 - uR - The field values at the right side of the face 2916 2917 Level: developer 2918 2919 .seealso: DMPlexGetFaceFields() 2920 @*/ 2921 PetscErrorCode DMPlexRestoreFaceFields(DM dm, PetscInt fStart, PetscInt fEnd, Vec locX, Vec locX_t, Vec faceGeometry, Vec cellGeometry, Vec locGrad, PetscInt *Nface, PetscScalar **uL, PetscScalar **uR) 2922 { 2923 PetscErrorCode ierr; 2924 2925 PetscFunctionBegin; 2926 ierr = DMRestoreWorkArray(dm, 0, MPIU_SCALAR, uL);CHKERRQ(ierr); 2927 ierr = DMRestoreWorkArray(dm, 0, MPIU_SCALAR, uR);CHKERRQ(ierr); 2928 PetscFunctionReturn(0); 2929 } 2930 2931 /*@C 2932 DMPlexGetFaceGeometry - Retrieve the geometric values for a chunk of faces 2933 2934 Input Parameters: 2935 + dm - The DM 2936 . fStart - The first face to include 2937 . fEnd - The first face to exclude 2938 . faceGeometry - A local vector with face geometry 2939 - cellGeometry - A local vector with cell geometry 2940 2941 Output Parameters: 2942 + Nface - The number of faces with field values 2943 . fgeom - The extract the face centroid and normal 2944 - vol - The cell volume 2945 2946 Level: developer 2947 2948 .seealso: DMPlexGetCellFields() 2949 @*/ 2950 PetscErrorCode DMPlexGetFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscInt *Nface, PetscFVFaceGeom **fgeom, PetscReal **vol) 2951 { 2952 DM dmFace, dmCell; 2953 DMLabel ghostLabel; 2954 const PetscScalar *facegeom, *cellgeom; 2955 PetscInt dim, numFaces = fEnd - fStart, iface, face; 2956 PetscErrorCode ierr; 2957 2958 PetscFunctionBegin; 2959 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2960 PetscValidHeaderSpecific(faceGeometry, VEC_CLASSID, 4); 2961 PetscValidHeaderSpecific(cellGeometry, VEC_CLASSID, 5); 2962 PetscValidPointer(fgeom, 6); 2963 PetscValidPointer(vol, 7); 2964 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2965 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 2966 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 2967 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 2968 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 2969 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 2970 ierr = PetscMalloc1(numFaces, fgeom);CHKERRQ(ierr); 2971 ierr = DMGetWorkArray(dm, numFaces*2, MPIU_SCALAR, vol);CHKERRQ(ierr); 2972 for (face = fStart, iface = 0; face < fEnd; ++face) { 2973 const PetscInt *cells; 2974 PetscFVFaceGeom *fg; 2975 PetscFVCellGeom *cgL, *cgR; 2976 PetscFVFaceGeom *fgeoml = *fgeom; 2977 PetscReal *voll = *vol; 2978 PetscInt ghost, d, nchild, nsupp; 2979 2980 ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr); 2981 ierr = DMPlexGetSupportSize(dm, face, &nsupp);CHKERRQ(ierr); 2982 ierr = DMPlexGetTreeChildren(dm, face, &nchild, NULL);CHKERRQ(ierr); 2983 if (ghost >= 0 || nsupp > 2 || nchild > 0) continue; 2984 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 2985 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 2986 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);CHKERRQ(ierr); 2987 ierr = DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);CHKERRQ(ierr); 2988 for (d = 0; d < dim; ++d) { 2989 fgeoml[iface].centroid[d] = fg->centroid[d]; 2990 fgeoml[iface].normal[d] = fg->normal[d]; 2991 } 2992 voll[iface*2+0] = cgL->volume; 2993 voll[iface*2+1] = cgR->volume; 2994 ++iface; 2995 } 2996 *Nface = iface; 2997 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 2998 ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 2999 PetscFunctionReturn(0); 3000 } 3001 3002 /*@C 3003 DMPlexRestoreFaceGeometry - Restore the field values values for a chunk of faces 3004 3005 Input Parameters: 3006 + dm - The DM 3007 . fStart - The first face to include 3008 . fEnd - The first face to exclude 3009 . faceGeometry - A local vector with face geometry 3010 - cellGeometry - A local vector with cell geometry 3011 3012 Output Parameters: 3013 + Nface - The number of faces with field values 3014 . fgeom - The extract the face centroid and normal 3015 - vol - The cell volume 3016 3017 Level: developer 3018 3019 .seealso: DMPlexGetFaceFields() 3020 @*/ 3021 PetscErrorCode DMPlexRestoreFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscInt *Nface, PetscFVFaceGeom **fgeom, PetscReal **vol) 3022 { 3023 PetscErrorCode ierr; 3024 3025 PetscFunctionBegin; 3026 ierr = PetscFree(*fgeom);CHKERRQ(ierr); 3027 ierr = DMRestoreWorkArray(dm, 0, MPIU_REAL, vol);CHKERRQ(ierr); 3028 PetscFunctionReturn(0); 3029 } 3030 3031 PetscErrorCode DMSNESGetFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom) 3032 { 3033 char composeStr[33] = {0}; 3034 PetscObjectId id; 3035 PetscContainer container; 3036 PetscErrorCode ierr; 3037 3038 PetscFunctionBegin; 3039 ierr = PetscObjectGetId((PetscObject)quad,&id);CHKERRQ(ierr); 3040 ierr = PetscSNPrintf(composeStr, 32, "DMSNESGetFEGeom_%x\n", id);CHKERRQ(ierr); 3041 ierr = PetscObjectQuery((PetscObject) pointIS, composeStr, (PetscObject *) &container);CHKERRQ(ierr); 3042 if (container) { 3043 ierr = PetscContainerGetPointer(container, (void **) geom);CHKERRQ(ierr); 3044 } else { 3045 ierr = DMFieldCreateFEGeom(coordField, pointIS, quad, faceData, geom);CHKERRQ(ierr); 3046 ierr = PetscContainerCreate(PETSC_COMM_SELF,&container);CHKERRQ(ierr); 3047 ierr = PetscContainerSetPointer(container, (void *) *geom);CHKERRQ(ierr); 3048 ierr = PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom);CHKERRQ(ierr); 3049 ierr = PetscObjectCompose((PetscObject) pointIS, composeStr, (PetscObject) container);CHKERRQ(ierr); 3050 ierr = PetscContainerDestroy(&container);CHKERRQ(ierr); 3051 } 3052 PetscFunctionReturn(0); 3053 } 3054 3055 PetscErrorCode DMSNESRestoreFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom) 3056 { 3057 PetscFunctionBegin; 3058 *geom = NULL; 3059 PetscFunctionReturn(0); 3060 } 3061 3062 PetscErrorCode DMPlexComputeResidual_Patch_Internal(DM dm, PetscSection section, IS cellIS, PetscReal t, Vec locX, Vec locX_t, Vec locF, void *user) 3063 { 3064 DM_Plex *mesh = (DM_Plex *) dm->data; 3065 const char *name = "Residual"; 3066 DM dmAux = NULL; 3067 DMLabel ghostLabel = NULL; 3068 PetscDS prob = NULL; 3069 PetscDS probAux = NULL; 3070 PetscBool useFEM = PETSC_FALSE; 3071 PetscBool isImplicit = (locX_t || t == PETSC_MIN_REAL) ? PETSC_TRUE : PETSC_FALSE; 3072 DMField coordField = NULL; 3073 Vec locA; 3074 PetscScalar *u = NULL, *u_t, *a, *uL = NULL, *uR = NULL; 3075 IS chunkIS; 3076 const PetscInt *cells; 3077 PetscInt cStart, cEnd, numCells; 3078 PetscInt Nf, f, totDim, totDimAux, numChunks, cellChunkSize, chunk, fStart, fEnd; 3079 PetscInt maxDegree = PETSC_MAX_INT; 3080 PetscQuadrature affineQuad = NULL, *quads = NULL; 3081 PetscFEGeom *affineGeom = NULL, **geoms = NULL; 3082 PetscErrorCode ierr; 3083 3084 PetscFunctionBegin; 3085 ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr); 3086 /* FEM+FVM */ 3087 /* 1: Get sizes from dm and dmAux */ 3088 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 3089 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 3090 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 3091 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 3092 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);CHKERRQ(ierr); 3093 if (locA) { 3094 ierr = VecGetDM(locA, &dmAux);CHKERRQ(ierr); 3095 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 3096 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 3097 } 3098 /* 2: Get geometric data */ 3099 for (f = 0; f < Nf; ++f) { 3100 PetscObject obj; 3101 PetscClassId id; 3102 PetscBool fimp; 3103 3104 ierr = PetscDSGetImplicit(prob, f, &fimp);CHKERRQ(ierr); 3105 if (isImplicit != fimp) continue; 3106 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 3107 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 3108 if (id == PETSCFE_CLASSID) {useFEM = PETSC_TRUE;} 3109 if (id == PETSCFV_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Use of FVM with PCPATCH not yet implemented"); 3110 } 3111 if (useFEM) { 3112 ierr = DMGetCoordinateField(dm, &coordField);CHKERRQ(ierr); 3113 ierr = DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);CHKERRQ(ierr); 3114 if (maxDegree <= 1) { 3115 ierr = DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);CHKERRQ(ierr); 3116 if (affineQuad) { 3117 ierr = DMSNESGetFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);CHKERRQ(ierr); 3118 } 3119 } else { 3120 ierr = PetscCalloc2(Nf,&quads,Nf,&geoms);CHKERRQ(ierr); 3121 for (f = 0; f < Nf; ++f) { 3122 PetscObject obj; 3123 PetscClassId id; 3124 PetscBool fimp; 3125 3126 ierr = PetscDSGetImplicit(prob, f, &fimp);CHKERRQ(ierr); 3127 if (isImplicit != fimp) continue; 3128 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 3129 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 3130 if (id == PETSCFE_CLASSID) { 3131 PetscFE fe = (PetscFE) obj; 3132 3133 ierr = PetscFEGetQuadrature(fe, &quads[f]);CHKERRQ(ierr); 3134 ierr = PetscObjectReference((PetscObject)quads[f]);CHKERRQ(ierr); 3135 ierr = DMSNESGetFEGeom(coordField,cellIS,quads[f],PETSC_FALSE,&geoms[f]);CHKERRQ(ierr); 3136 } 3137 } 3138 } 3139 } 3140 /* Loop over chunks */ 3141 ierr = ISGetPointRange(cellIS, &cStart, &cEnd, &cells);CHKERRQ(ierr); 3142 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 3143 if (useFEM) {ierr = ISCreate(PETSC_COMM_SELF, &chunkIS);CHKERRQ(ierr);} 3144 numCells = cEnd - cStart; 3145 numChunks = 1; 3146 cellChunkSize = numCells/numChunks; 3147 numChunks = PetscMin(1,numCells); 3148 for (chunk = 0; chunk < numChunks; ++chunk) { 3149 PetscScalar *elemVec, *fluxL = NULL, *fluxR = NULL; 3150 PetscReal *vol = NULL; 3151 PetscFVFaceGeom *fgeom = NULL; 3152 PetscInt cS = cStart+chunk*cellChunkSize, cE = PetscMin(cS+cellChunkSize, cEnd), numCells = cE - cS, c; 3153 PetscInt numFaces = 0; 3154 3155 /* Extract field coefficients */ 3156 if (useFEM) { 3157 ierr = ISGetPointSubrange(chunkIS, cS, cE, cells);CHKERRQ(ierr); 3158 ierr = DMPlexGetCellFields(dm, chunkIS, locX, locX_t, locA, &u, &u_t, &a);CHKERRQ(ierr); 3159 ierr = DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);CHKERRQ(ierr); 3160 ierr = PetscMemzero(elemVec, numCells*totDim * sizeof(PetscScalar));CHKERRQ(ierr); 3161 } 3162 /* TODO We will interlace both our field coefficients (u, u_t, uL, uR, etc.) and our output (elemVec, fL, fR). I think this works */ 3163 /* Loop over fields */ 3164 for (f = 0; f < Nf; ++f) { 3165 PetscObject obj; 3166 PetscClassId id; 3167 PetscBool fimp; 3168 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 3169 3170 ierr = PetscDSGetImplicit(prob, f, &fimp);CHKERRQ(ierr); 3171 if (isImplicit != fimp) continue; 3172 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 3173 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 3174 if (id == PETSCFE_CLASSID) { 3175 PetscFE fe = (PetscFE) obj; 3176 PetscFEGeom *geom = affineGeom ? affineGeom : geoms[f]; 3177 PetscFEGeom *chunkGeom = NULL; 3178 PetscQuadrature quad = affineQuad ? affineQuad : quads[f]; 3179 PetscInt Nq, Nb; 3180 3181 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 3182 ierr = PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 3183 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 3184 blockSize = Nb; 3185 batchSize = numBlocks * blockSize; 3186 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 3187 numChunks = numCells / (numBatches*batchSize); 3188 Ne = numChunks*numBatches*batchSize; 3189 Nr = numCells % (numBatches*batchSize); 3190 offset = numCells - Nr; 3191 /* Integrate FE residual to get elemVec (need fields at quadrature points) */ 3192 /* For FV, I think we use a P0 basis and the cell coefficients (for subdivided cells, we can tweak the basis tabulation to be the indicator function) */ 3193 ierr = PetscFEGeomGetChunk(geom,0,offset,&chunkGeom);CHKERRQ(ierr); 3194 ierr = PetscFEIntegrateResidual(fe, prob, f, Ne, chunkGeom, u, u_t, probAux, a, t, elemVec);CHKERRQ(ierr); 3195 ierr = PetscFEGeomGetChunk(geom,offset,numCells,&chunkGeom);CHKERRQ(ierr); 3196 ierr = PetscFEIntegrateResidual(fe, prob, f, Nr, chunkGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, &elemVec[offset*totDim]);CHKERRQ(ierr); 3197 ierr = PetscFEGeomRestoreChunk(geom,offset,numCells,&chunkGeom);CHKERRQ(ierr); 3198 } else if (id == PETSCFV_CLASSID) { 3199 PetscFV fv = (PetscFV) obj; 3200 3201 Ne = numFaces; 3202 /* Riemann solve over faces (need fields at face centroids) */ 3203 /* We need to evaluate FE fields at those coordinates */ 3204 ierr = PetscFVIntegrateRHSFunction(fv, prob, f, Ne, fgeom, vol, uL, uR, fluxL, fluxR);CHKERRQ(ierr); 3205 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 3206 } 3207 /* Loop over domain */ 3208 if (useFEM) { 3209 /* Add elemVec to locX */ 3210 for (c = cS; c < cE; ++c) { 3211 const PetscInt cell = cells ? cells[c] : c; 3212 const PetscInt cind = c - cStart; 3213 3214 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(cell, name, totDim, &elemVec[cind*totDim]);CHKERRQ(ierr);} 3215 if (ghostLabel) { 3216 PetscInt ghostVal; 3217 3218 ierr = DMLabelGetValue(ghostLabel,cell,&ghostVal);CHKERRQ(ierr); 3219 if (ghostVal > 0) continue; 3220 } 3221 ierr = DMPlexVecSetClosure(dm, section, locF, cell, &elemVec[cind*totDim], ADD_ALL_VALUES);CHKERRQ(ierr); 3222 } 3223 } 3224 /* Handle time derivative */ 3225 if (locX_t) { 3226 PetscScalar *x_t, *fa; 3227 3228 ierr = VecGetArray(locF, &fa);CHKERRQ(ierr); 3229 ierr = VecGetArray(locX_t, &x_t);CHKERRQ(ierr); 3230 for (f = 0; f < Nf; ++f) { 3231 PetscFV fv; 3232 PetscObject obj; 3233 PetscClassId id; 3234 PetscInt pdim, d; 3235 3236 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 3237 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 3238 if (id != PETSCFV_CLASSID) continue; 3239 fv = (PetscFV) obj; 3240 ierr = PetscFVGetNumComponents(fv, &pdim);CHKERRQ(ierr); 3241 for (c = cS; c < cE; ++c) { 3242 const PetscInt cell = cells ? cells[c] : c; 3243 PetscScalar *u_t, *r; 3244 3245 if (ghostLabel) { 3246 PetscInt ghostVal; 3247 3248 ierr = DMLabelGetValue(ghostLabel, cell, &ghostVal);CHKERRQ(ierr); 3249 if (ghostVal > 0) continue; 3250 } 3251 ierr = DMPlexPointLocalFieldRead(dm, cell, f, x_t, &u_t);CHKERRQ(ierr); 3252 ierr = DMPlexPointLocalFieldRef(dm, cell, f, fa, &r);CHKERRQ(ierr); 3253 for (d = 0; d < pdim; ++d) r[d] += u_t[d]; 3254 } 3255 } 3256 ierr = VecRestoreArray(locX_t, &x_t);CHKERRQ(ierr); 3257 ierr = VecRestoreArray(locF, &fa);CHKERRQ(ierr); 3258 } 3259 if (useFEM) { 3260 ierr = DMPlexRestoreCellFields(dm, chunkIS, locX, locX_t, locA, &u, &u_t, &a);CHKERRQ(ierr); 3261 ierr = DMRestoreWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);CHKERRQ(ierr); 3262 } 3263 } 3264 if (useFEM) {ierr = ISDestroy(&chunkIS);CHKERRQ(ierr);} 3265 ierr = ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);CHKERRQ(ierr); 3266 /* TODO Could include boundary residual here (see DMPlexComputeResidual_Internal) */ 3267 if (useFEM) { 3268 if (maxDegree <= 1) { 3269 ierr = DMSNESRestoreFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);CHKERRQ(ierr); 3270 ierr = PetscQuadratureDestroy(&affineQuad);CHKERRQ(ierr); 3271 } else { 3272 for (f = 0; f < Nf; ++f) { 3273 ierr = DMSNESRestoreFEGeom(coordField,cellIS,quads[f],PETSC_FALSE,&geoms[f]);CHKERRQ(ierr); 3274 ierr = PetscQuadratureDestroy(&quads[f]);CHKERRQ(ierr); 3275 } 3276 ierr = PetscFree2(quads,geoms);CHKERRQ(ierr); 3277 } 3278 } 3279 ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr); 3280 PetscFunctionReturn(0); 3281 } 3282 3283 /* 3284 We always assemble JacP, and if the matrix is different from Jac and two different sets of point functions are provided, we also assemble Jac 3285 3286 X - The local solution vector 3287 X_t - The local solution time derviative vector, or NULL 3288 */ 3289 PetscErrorCode DMPlexComputeJacobian_Patch_Internal(DM dm, PetscSection section, PetscSection globalSection, IS cellIS, 3290 PetscReal t, PetscReal X_tShift, Vec X, Vec X_t, Mat Jac, Mat JacP, void *ctx) 3291 { 3292 DM_Plex *mesh = (DM_Plex *) dm->data; 3293 const char *name = "Jacobian", *nameP = "JacobianPre"; 3294 DM dmAux = NULL; 3295 PetscDS prob, probAux = NULL; 3296 PetscSection sectionAux = NULL; 3297 Vec A; 3298 DMField coordField; 3299 PetscFEGeom *cgeomFEM; 3300 PetscQuadrature qGeom = NULL; 3301 Mat J = Jac, JP = JacP; 3302 PetscScalar *work, *u = NULL, *u_t = NULL, *a = NULL, *elemMat = NULL, *elemMatP = NULL, *elemMatD = NULL; 3303 PetscBool hasJac, hasPrec, hasDyn, assembleJac, isMatIS, isMatISP, *isFE, hasFV = PETSC_FALSE; 3304 const PetscInt *cells; 3305 PetscInt Nf, fieldI, fieldJ, maxDegree, numCells, cStart, cEnd, numChunks, chunkSize, chunk, totDim, totDimAux = 0, sz, wsz, off = 0, offCell = 0; 3306 PetscErrorCode ierr; 3307 3308 PetscFunctionBegin; 3309 CHKMEMQ; 3310 ierr = ISGetLocalSize(cellIS, &numCells);CHKERRQ(ierr); 3311 ierr = ISGetPointRange(cellIS, &cStart, &cEnd, &cells);CHKERRQ(ierr); 3312 ierr = PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr); 3313 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 3314 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 3315 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 3316 if (dmAux) { 3317 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 3318 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 3319 } 3320 /* Get flags */ 3321 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 3322 ierr = DMGetWorkArray(dm, Nf, MPIU_BOOL, &isFE);CHKERRQ(ierr); 3323 for (fieldI = 0; fieldI < Nf; ++fieldI) { 3324 PetscObject disc; 3325 PetscClassId id; 3326 ierr = PetscDSGetDiscretization(prob, fieldI, &disc);CHKERRQ(ierr); 3327 ierr = PetscObjectGetClassId(disc, &id);CHKERRQ(ierr); 3328 if (id == PETSCFE_CLASSID) {isFE[fieldI] = PETSC_TRUE;} 3329 else if (id == PETSCFV_CLASSID) {hasFV = PETSC_TRUE; isFE[fieldI] = PETSC_FALSE;} 3330 } 3331 ierr = PetscDSHasJacobian(prob, &hasJac);CHKERRQ(ierr); 3332 ierr = PetscDSHasJacobianPreconditioner(prob, &hasPrec);CHKERRQ(ierr); 3333 ierr = PetscDSHasDynamicJacobian(prob, &hasDyn);CHKERRQ(ierr); 3334 assembleJac = hasJac && hasPrec && (Jac != JacP) ? PETSC_TRUE : PETSC_FALSE; 3335 hasDyn = hasDyn && (X_tShift != 0.0) ? PETSC_TRUE : PETSC_FALSE; 3336 ierr = PetscObjectTypeCompare((PetscObject) Jac, MATIS, &isMatIS);CHKERRQ(ierr); 3337 ierr = PetscObjectTypeCompare((PetscObject) JacP, MATIS, &isMatISP);CHKERRQ(ierr); 3338 /* Setup input data and temp arrays (should be DMGetWorkArray) */ 3339 if (isMatISP || isMatISP) {ierr = DMPlexGetSubdomainSection(dm, &globalSection);CHKERRQ(ierr);} 3340 if (isMatIS) {ierr = MatISGetLocalMat(Jac, &J);CHKERRQ(ierr);} 3341 if (isMatISP) {ierr = MatISGetLocalMat(JacP, &JP);CHKERRQ(ierr);} 3342 if (hasFV) {ierr = MatSetOption(JP, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);CHKERRQ(ierr);} /* No allocated space for FV stuff, so ignore the zero entries */ 3343 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 3344 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 3345 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 3346 if (probAux) {ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr);} 3347 CHKMEMQ; 3348 /* Compute batch sizes */ 3349 if (isFE[0]) { 3350 PetscFE fe; 3351 PetscQuadrature q; 3352 PetscInt numQuadPoints, numBatches, batchSize, numBlocks, blockSize, Nb; 3353 3354 ierr = PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);CHKERRQ(ierr); 3355 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 3356 ierr = PetscQuadratureGetData(q, NULL, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 3357 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 3358 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 3359 blockSize = Nb*numQuadPoints; 3360 batchSize = numBlocks * blockSize; 3361 chunkSize = numBatches * batchSize; 3362 numChunks = numCells / chunkSize + numCells % chunkSize; 3363 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 3364 } else { 3365 chunkSize = numCells; 3366 numChunks = 1; 3367 } 3368 /* Get work space */ 3369 wsz = (((X?1:0) + (X_t?1:0) + (dmAux?1:0))*totDim + ((hasJac?1:0) + (hasPrec?1:0) + (hasDyn?1:0))*totDim*totDim)*chunkSize; 3370 ierr = DMGetWorkArray(dm, wsz, MPIU_SCALAR, &work);CHKERRQ(ierr); 3371 ierr = PetscMemzero(work, wsz * sizeof(PetscScalar));CHKERRQ(ierr); 3372 off = 0; 3373 u = X ? (sz = chunkSize*totDim, off += sz, work+off-sz) : NULL; 3374 u_t = X_t ? (sz = chunkSize*totDim, off += sz, work+off-sz) : NULL; 3375 a = dmAux ? (sz = chunkSize*totDimAux, off += sz, work+off-sz) : NULL; 3376 elemMat = hasJac ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL; 3377 elemMatP = hasPrec ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL; 3378 elemMatD = hasDyn ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL; 3379 if (off != wsz) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error is workspace size %D should be %D", off, wsz); 3380 /* Setup geometry */ 3381 ierr = DMGetCoordinateField(dm, &coordField);CHKERRQ(ierr); 3382 ierr = DMFieldGetDegree(coordField, cellIS, NULL, &maxDegree);CHKERRQ(ierr); 3383 if (maxDegree <= 1) {ierr = DMFieldCreateDefaultQuadrature(coordField, cellIS, &qGeom);CHKERRQ(ierr);} 3384 if (!qGeom) { 3385 PetscFE fe; 3386 3387 ierr = PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);CHKERRQ(ierr); 3388 ierr = PetscFEGetQuadrature(fe, &qGeom);CHKERRQ(ierr); 3389 ierr = PetscObjectReference((PetscObject) qGeom);CHKERRQ(ierr); 3390 } 3391 ierr = DMSNESGetFEGeom(coordField, cellIS, qGeom, PETSC_FALSE, &cgeomFEM);CHKERRQ(ierr); 3392 /* Compute volume integrals */ 3393 if (assembleJac) {ierr = MatZeroEntries(J);CHKERRQ(ierr);} 3394 ierr = MatZeroEntries(JP);CHKERRQ(ierr); 3395 for (chunk = 0; chunk < numChunks; ++chunk, offCell += chunkSize) { 3396 const PetscInt Ncell = PetscMin(chunkSize, numCells - offCell); 3397 PetscInt c; 3398 3399 /* Extract values */ 3400 for (c = 0; c < Ncell; ++c) { 3401 const PetscInt cell = cells ? cells[c+offCell] : c+offCell; 3402 PetscScalar *x = NULL, *x_t = NULL; 3403 PetscInt i; 3404 3405 if (X) { 3406 ierr = DMPlexVecGetClosure(dm, section, X, cell, NULL, &x);CHKERRQ(ierr); 3407 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 3408 ierr = DMPlexVecRestoreClosure(dm, section, X, cell, NULL, &x);CHKERRQ(ierr); 3409 } 3410 if (X_t) { 3411 ierr = DMPlexVecGetClosure(dm, section, X_t, cell, NULL, &x_t);CHKERRQ(ierr); 3412 for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i]; 3413 ierr = DMPlexVecRestoreClosure(dm, section, X_t, cell, NULL, &x_t);CHKERRQ(ierr); 3414 } 3415 if (dmAux) { 3416 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, cell, NULL, &x);CHKERRQ(ierr); 3417 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 3418 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, cell, NULL, &x);CHKERRQ(ierr); 3419 } 3420 } 3421 CHKMEMQ; 3422 for (fieldI = 0; fieldI < Nf; ++fieldI) { 3423 PetscFE fe; 3424 ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr); 3425 for (fieldJ = 0; fieldJ < Nf; ++fieldJ) { 3426 if (hasJac) {ierr = PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMat);CHKERRQ(ierr);} 3427 if (hasPrec) {ierr = PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_PRE, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMatP);CHKERRQ(ierr);} 3428 if (hasDyn) {ierr = PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMatD);CHKERRQ(ierr);} 3429 } 3430 /* For finite volume, add the identity */ 3431 if (!isFE[fieldI]) { 3432 PetscFV fv; 3433 PetscInt eOffset = 0, Nc, fc, foff; 3434 3435 ierr = PetscDSGetFieldOffset(prob, fieldI, &foff);CHKERRQ(ierr); 3436 ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fv);CHKERRQ(ierr); 3437 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 3438 for (c = 0; c < chunkSize; ++c, eOffset += totDim*totDim) { 3439 for (fc = 0; fc < Nc; ++fc) { 3440 const PetscInt i = foff + fc; 3441 if (hasJac) {elemMat [eOffset+i*totDim+i] = 1.0;} 3442 if (hasPrec) {elemMatP[eOffset+i*totDim+i] = 1.0;} 3443 } 3444 } 3445 } 3446 } 3447 CHKMEMQ; 3448 /* Add contribution from X_t */ 3449 if (hasDyn) {for (c = 0; c < chunkSize*totDim*totDim; ++c) elemMat[c] += X_tShift*elemMatD[c];} 3450 /* Insert values into matrix */ 3451 for (c = 0; c < Ncell; ++c) { 3452 const PetscInt cell = cells ? cells[c+offCell] : c+offCell; 3453 if (mesh->printFEM > 1) { 3454 if (hasJac) {ierr = DMPrintCellMatrix(cell, name, totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);CHKERRQ(ierr);} 3455 if (hasPrec) {ierr = DMPrintCellMatrix(cell, nameP, totDim, totDim, &elemMatP[(c-cStart)*totDim*totDim]);CHKERRQ(ierr);} 3456 } 3457 if (assembleJac) {ierr = DMPlexMatSetClosure(dm, section, globalSection, Jac, cell, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);CHKERRQ(ierr);} 3458 ierr = DMPlexMatSetClosure(dm, section, globalSection, JP, cell, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);CHKERRQ(ierr); 3459 } 3460 CHKMEMQ; 3461 } 3462 /* Cleanup */ 3463 ierr = DMSNESRestoreFEGeom(coordField, cellIS, qGeom, PETSC_FALSE, &cgeomFEM);CHKERRQ(ierr); 3464 ierr = PetscQuadratureDestroy(&qGeom);CHKERRQ(ierr); 3465 if (hasFV) {ierr = MatSetOption(JacP, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);CHKERRQ(ierr);} 3466 ierr = DMRestoreWorkArray(dm, Nf, MPIU_BOOL, &isFE);CHKERRQ(ierr); 3467 ierr = DMRestoreWorkArray(dm, ((1 + (X_t?1:0) + (dmAux?1:0))*totDim + ((hasJac?1:0) + (hasPrec?1:0) + (hasDyn?1:0))*totDim*totDim)*chunkSize, MPIU_SCALAR, &work);CHKERRQ(ierr); 3468 /* Compute boundary integrals */ 3469 /* ierr = DMPlexComputeBdJacobian_Internal(dm, X, X_t, t, X_tShift, Jac, JacP, ctx);CHKERRQ(ierr); */ 3470 /* Assemble matrix */ 3471 if (assembleJac) {ierr = MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);ierr = MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);} 3472 ierr = MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);ierr = MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3473 ierr = PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr); 3474 CHKMEMQ; 3475 PetscFunctionReturn(0); 3476 } 3477