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