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