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