1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petscsf.h> 3 4 #include <petsc/private/petscfeimpl.h> 5 #include <petsc/private/petscfvimpl.h> 6 7 /*@ 8 DMPlexGetScale - Get the scale for the specified fundamental unit 9 10 Not collective 11 12 Input Arguments: 13 + dm - the DM 14 - unit - The SI unit 15 16 Output Argument: 17 . scale - The value used to scale all quantities with this unit 18 19 Level: advanced 20 21 .seealso: DMPlexSetScale(), PetscUnit 22 @*/ 23 PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale) 24 { 25 DM_Plex *mesh = (DM_Plex*) dm->data; 26 27 PetscFunctionBegin; 28 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 29 PetscValidPointer(scale, 3); 30 *scale = mesh->scale[unit]; 31 PetscFunctionReturn(0); 32 } 33 34 /*@ 35 DMPlexSetScale - Set the scale for the specified fundamental unit 36 37 Not collective 38 39 Input Arguments: 40 + dm - the DM 41 . unit - The SI unit 42 - scale - The value used to scale all quantities with this unit 43 44 Level: advanced 45 46 .seealso: DMPlexGetScale(), PetscUnit 47 @*/ 48 PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale) 49 { 50 DM_Plex *mesh = (DM_Plex*) dm->data; 51 52 PetscFunctionBegin; 53 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 54 mesh->scale[unit] = scale; 55 PetscFunctionReturn(0); 56 } 57 58 static PetscErrorCode DMPlexProjectRigidBody_Private(PetscInt dim, PetscReal t, const PetscReal X[], PetscInt Nf, PetscScalar *mode, void *ctx) 59 { 60 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}}}; 61 PetscInt *ctxInt = (PetscInt *) ctx; 62 PetscInt dim2 = ctxInt[0]; 63 PetscInt d = ctxInt[1]; 64 PetscInt i, j, k = dim > 2 ? d - dim : d; 65 66 PetscFunctionBegin; 67 if (dim != dim2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input dimension %d does not match context dimension %d", dim, dim2); 68 for (i = 0; i < dim; i++) mode[i] = 0.; 69 if (d < dim) { 70 mode[d] = 1.; /* Translation along axis d */ 71 } else { 72 for (i = 0; i < dim; i++) { 73 for (j = 0; j < dim; j++) { 74 mode[j] += eps[i][j][k]*X[i]; /* Rotation about axis d */ 75 } 76 } 77 } 78 PetscFunctionReturn(0); 79 } 80 81 /*@C 82 DMPlexCreateRigidBody - For the default global section, create rigid body modes by function space interpolation 83 84 Collective on DM 85 86 Input Arguments: 87 . dm - the DM 88 89 Output Argument: 90 . sp - the null space 91 92 Note: This is necessary to provide a suitable coarse space for algebraic multigrid 93 94 Level: advanced 95 96 .seealso: MatNullSpaceCreate(), PCGAMG 97 @*/ 98 PetscErrorCode DMPlexCreateRigidBody(DM dm, MatNullSpace *sp) 99 { 100 MPI_Comm comm; 101 Vec mode[6]; 102 PetscSection section, globalSection; 103 PetscInt dim, dimEmbed, n, m, d, i, j; 104 PetscErrorCode ierr; 105 106 PetscFunctionBegin; 107 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 108 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 109 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 110 if (dim == 1) { 111 ierr = MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);CHKERRQ(ierr); 112 PetscFunctionReturn(0); 113 } 114 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 115 ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr); 116 ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr); 117 m = (dim*(dim+1))/2; 118 ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr); 119 ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr); 120 ierr = VecSetUp(mode[0]);CHKERRQ(ierr); 121 for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);} 122 for (d = 0; d < m; d++) { 123 PetscInt ctx[2]; 124 PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private; 125 void *voidctx = (void *) (&ctx[0]); 126 127 ctx[0] = dimEmbed; 128 ctx[1] = d; 129 ierr = DMProjectFunction(dm, 0.0, &func, &voidctx, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 130 } 131 for (i = 0; i < dim; ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);} 132 /* Orthonormalize system */ 133 for (i = dim; i < m; ++i) { 134 PetscScalar dots[6]; 135 136 ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr); 137 for (j = 0; j < i; ++j) dots[j] *= -1.0; 138 ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr); 139 ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr); 140 } 141 ierr = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);CHKERRQ(ierr); 142 for (i = 0; i< m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);} 143 PetscFunctionReturn(0); 144 } 145 146 /*@C 147 DMPlexCreateRigidBodies - For the default global section, create rigid body modes by function space interpolation 148 149 Collective on DM 150 151 Input Arguments: 152 + dm - the DM 153 . nb - The number of bodies 154 . label - The DMLabel marking each domain 155 . nids - The number of ids per body 156 - ids - An array of the label ids in sequence for each domain 157 158 Output Argument: 159 . sp - the null space 160 161 Note: This is necessary to provide a suitable coarse space for algebraic multigrid 162 163 Level: advanced 164 165 .seealso: MatNullSpaceCreate() 166 @*/ 167 PetscErrorCode DMPlexCreateRigidBodies(DM dm, PetscInt nb, DMLabel label, const PetscInt nids[], const PetscInt ids[], MatNullSpace *sp) 168 { 169 MPI_Comm comm; 170 PetscSection section, globalSection; 171 Vec *mode; 172 PetscScalar *dots; 173 PetscInt dim, dimEmbed, n, m, b, d, i, j, off; 174 PetscErrorCode ierr; 175 176 PetscFunctionBegin; 177 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 178 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 179 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 180 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 181 ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr); 182 ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr); 183 m = nb * (dim*(dim+1))/2; 184 ierr = PetscMalloc2(m, &mode, m, &dots);CHKERRQ(ierr); 185 ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr); 186 ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr); 187 ierr = VecSetUp(mode[0]);CHKERRQ(ierr); 188 for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);} 189 for (b = 0, off = 0; b < nb; ++b) { 190 for (d = 0; d < m/nb; ++d) { 191 PetscInt ctx[2]; 192 PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private; 193 void *voidctx = (void *) (&ctx[0]); 194 195 ctx[0] = dimEmbed; 196 ctx[1] = d; 197 ierr = DMProjectFunctionLabel(dm, 0.0, label, nids[b], &ids[off], 0, NULL, &func, &voidctx, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 198 off += nids[b]; 199 } 200 } 201 for (i = 0; i < dim; ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);} 202 /* Orthonormalize system */ 203 for (i = 0; i < m; ++i) { 204 ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr); 205 for (j = 0; j < i; ++j) dots[j] *= -1.0; 206 ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr); 207 ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr); 208 } 209 ierr = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);CHKERRQ(ierr); 210 for (i = 0; i< m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);} 211 ierr = PetscFree2(mode, dots);CHKERRQ(ierr); 212 PetscFunctionReturn(0); 213 } 214 215 /*@ 216 DMPlexSetMaxProjectionHeight - In DMPlexProjectXXXLocal() functions, the projected values of a basis function's dofs 217 are computed by associating the basis function with one of the mesh points in its transitively-closed support, and 218 evaluating the dual space basis of that point. A basis function is associated with the point in its 219 transitively-closed support whose mesh height is highest (w.r.t. DAG height), but not greater than the maximum 220 projection height, which is set with this function. By default, the maximum projection height is zero, which means 221 that only mesh cells are used to project basis functions. A height of one, for example, evaluates a cell-interior 222 basis functions using its cells dual space basis, but all other basis functions with the dual space basis of a face. 223 224 Input Parameters: 225 + dm - the DMPlex object 226 - height - the maximum projection height >= 0 227 228 Level: advanced 229 230 .seealso: DMPlexGetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal() 231 @*/ 232 PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height) 233 { 234 DM_Plex *plex = (DM_Plex *) dm->data; 235 236 PetscFunctionBegin; 237 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 238 plex->maxProjectionHeight = height; 239 PetscFunctionReturn(0); 240 } 241 242 /*@ 243 DMPlexGetMaxProjectionHeight - Get the maximum height (w.r.t. DAG) of mesh points used to evaluate dual bases in 244 DMPlexProjectXXXLocal() functions. 245 246 Input Parameters: 247 . dm - the DMPlex object 248 249 Output Parameters: 250 . height - the maximum projection height 251 252 Level: intermediate 253 254 .seealso: DMPlexSetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal() 255 @*/ 256 PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height) 257 { 258 DM_Plex *plex = (DM_Plex *) dm->data; 259 260 PetscFunctionBegin; 261 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 262 *height = plex->maxProjectionHeight; 263 PetscFunctionReturn(0); 264 } 265 266 /*@C 267 DMPlexInsertBoundaryValuesEssential - Insert boundary values into a local vector 268 269 Input Parameters: 270 + dm - The DM, with a PetscDS that matches the problem being constrained 271 . time - The time 272 . field - The field to constrain 273 . Nc - The number of constrained field components, or 0 for all components 274 . comps - An array of constrained component numbers, or NULL for all components 275 . label - The DMLabel defining constrained points 276 . numids - The number of DMLabel ids for constrained points 277 . ids - An array of ids for constrained points 278 . func - A pointwise function giving boundary values 279 - ctx - An optional user context for bcFunc 280 281 Output Parameter: 282 . locX - A local vector to receives the boundary values 283 284 Level: developer 285 286 .seealso: DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary() 287 @*/ 288 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) 289 { 290 PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx); 291 void **ctxs; 292 PetscInt numFields; 293 PetscErrorCode ierr; 294 295 PetscFunctionBegin; 296 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 297 ierr = PetscCalloc2(numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 298 funcs[field] = func; 299 ctxs[field] = ctx; 300 ierr = DMProjectFunctionLabelLocal(dm, time, label, numids, ids, Nc, comps, funcs, ctxs, INSERT_BC_VALUES, locX);CHKERRQ(ierr); 301 ierr = PetscFree2(funcs,ctxs);CHKERRQ(ierr); 302 PetscFunctionReturn(0); 303 } 304 305 /*@C 306 DMPlexInsertBoundaryValuesEssentialField - Insert boundary values into a local vector 307 308 Input Parameters: 309 + dm - The DM, with a PetscDS that matches the problem being constrained 310 . time - The time 311 . locU - A local vector with the input solution values 312 . field - The field to constrain 313 . Nc - The number of constrained field components, or 0 for all components 314 . comps - An array of constrained component numbers, or NULL for all components 315 . label - The DMLabel defining constrained points 316 . numids - The number of DMLabel ids for constrained points 317 . ids - An array of ids for constrained points 318 . func - A pointwise function giving boundary values 319 - ctx - An optional user context for bcFunc 320 321 Output Parameter: 322 . locX - A local vector to receives the boundary values 323 324 Level: developer 325 326 .seealso: DMPlexInsertBoundaryValuesEssential(), DMAddBoundary() 327 @*/ 328 PetscErrorCode DMPlexInsertBoundaryValuesEssentialField(DM dm, PetscReal time, Vec locU, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[], 329 void (*func)(PetscInt, PetscInt, PetscInt, 330 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 331 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 332 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], 333 PetscScalar[]), 334 void *ctx, Vec locX) 335 { 336 void (**funcs)(PetscInt, PetscInt, PetscInt, 337 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 338 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 339 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]); 340 void **ctxs; 341 PetscInt numFields; 342 PetscErrorCode ierr; 343 344 PetscFunctionBegin; 345 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 346 ierr = PetscCalloc2(numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 347 funcs[field] = func; 348 ctxs[field] = ctx; 349 ierr = DMProjectFieldLabelLocal(dm, time, label, numids, ids, Nc, comps, locU, funcs, INSERT_BC_VALUES, locX);CHKERRQ(ierr); 350 ierr = PetscFree2(funcs,ctxs);CHKERRQ(ierr); 351 PetscFunctionReturn(0); 352 } 353 354 /*@C 355 DMPlexInsertBoundaryValuesRiemann - Insert boundary values into a local vector 356 357 Input Parameters: 358 + dm - The DM, with a PetscDS that matches the problem being constrained 359 . time - The time 360 . faceGeometry - A vector with the FVM face geometry information 361 . cellGeometry - A vector with the FVM cell geometry information 362 . Grad - A vector with the FVM cell gradient information 363 . field - The field to constrain 364 . Nc - The number of constrained field components, or 0 for all components 365 . comps - An array of constrained component numbers, or NULL for all components 366 . label - The DMLabel defining constrained points 367 . numids - The number of DMLabel ids for constrained points 368 . ids - An array of ids for constrained points 369 . func - A pointwise function giving boundary values 370 - ctx - An optional user context for bcFunc 371 372 Output Parameter: 373 . locX - A local vector to receives the boundary values 374 375 Note: This implementation currently ignores the numcomps/comps argument from DMAddBoundary() 376 377 Level: developer 378 379 .seealso: DMPlexInsertBoundaryValuesEssential(), DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary() 380 @*/ 381 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[], 382 PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX) 383 { 384 PetscDS prob; 385 PetscSF sf; 386 DM dmFace, dmCell, dmGrad; 387 const PetscScalar *facegeom, *cellgeom = NULL, *grad; 388 const PetscInt *leaves; 389 PetscScalar *x, *fx; 390 PetscInt dim, nleaves, loc, fStart, fEnd, pdim, i; 391 PetscErrorCode ierr, ierru = 0; 392 393 PetscFunctionBegin; 394 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 395 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 396 nleaves = PetscMax(0, nleaves); 397 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 398 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 399 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 400 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 401 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 402 if (cellGeometry) { 403 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 404 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 405 } 406 if (Grad) { 407 PetscFV fv; 408 409 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);CHKERRQ(ierr); 410 ierr = VecGetDM(Grad, &dmGrad);CHKERRQ(ierr); 411 ierr = VecGetArrayRead(Grad, &grad);CHKERRQ(ierr); 412 ierr = PetscFVGetNumComponents(fv, &pdim);CHKERRQ(ierr); 413 ierr = DMGetWorkArray(dm, pdim, MPIU_SCALAR, &fx);CHKERRQ(ierr); 414 } 415 ierr = VecGetArray(locX, &x);CHKERRQ(ierr); 416 for (i = 0; i < numids; ++i) { 417 IS faceIS; 418 const PetscInt *faces; 419 PetscInt numFaces, f; 420 421 ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr); 422 if (!faceIS) continue; /* No points with that id on this process */ 423 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 424 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 425 for (f = 0; f < numFaces; ++f) { 426 const PetscInt face = faces[f], *cells; 427 PetscFVFaceGeom *fg; 428 429 if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */ 430 ierr = PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);CHKERRQ(ierr); 431 if (loc >= 0) continue; 432 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 433 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 434 if (Grad) { 435 PetscFVCellGeom *cg; 436 PetscScalar *cx, *cgrad; 437 PetscScalar *xG; 438 PetscReal dx[3]; 439 PetscInt d; 440 441 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);CHKERRQ(ierr); 442 ierr = DMPlexPointLocalRead(dm, cells[0], x, &cx);CHKERRQ(ierr); 443 ierr = DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);CHKERRQ(ierr); 444 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 445 DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx); 446 for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx); 447 ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx); 448 if (ierru) { 449 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 450 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 451 goto cleanup; 452 } 453 } else { 454 PetscScalar *xI; 455 PetscScalar *xG; 456 457 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr); 458 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 459 ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx); 460 if (ierru) { 461 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 462 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 463 goto cleanup; 464 } 465 } 466 } 467 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 468 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 469 } 470 cleanup: 471 ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr); 472 if (Grad) { 473 ierr = DMRestoreWorkArray(dm, pdim, MPIU_SCALAR, &fx);CHKERRQ(ierr); 474 ierr = VecRestoreArrayRead(Grad, &grad);CHKERRQ(ierr); 475 } 476 if (cellGeometry) {ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr);} 477 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 478 CHKERRQ(ierru); 479 PetscFunctionReturn(0); 480 } 481 482 PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM) 483 { 484 PetscInt numBd, b; 485 PetscErrorCode ierr; 486 487 PetscFunctionBegin; 488 ierr = PetscDSGetNumBoundary(dm->prob, &numBd);CHKERRQ(ierr); 489 for (b = 0; b < numBd; ++b) { 490 DMBoundaryConditionType type; 491 const char *labelname; 492 DMLabel label; 493 PetscInt field, Nc; 494 const PetscInt *comps; 495 PetscObject obj; 496 PetscClassId id; 497 void (*func)(void); 498 PetscInt numids; 499 const PetscInt *ids; 500 void *ctx; 501 502 ierr = DMGetBoundary(dm, b, &type, NULL, &labelname, &field, &Nc, &comps, &func, &numids, &ids, &ctx);CHKERRQ(ierr); 503 if (insertEssential != (type & DM_BC_ESSENTIAL)) continue; 504 ierr = DMGetLabel(dm, labelname, &label);CHKERRQ(ierr); 505 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 506 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 507 if (id == PETSCFE_CLASSID) { 508 switch (type) { 509 /* for FEM, there is no insertion to be done for non-essential boundary conditions */ 510 case DM_BC_ESSENTIAL: 511 ierr = DMPlexLabelAddCells(dm,label);CHKERRQ(ierr); 512 ierr = DMPlexInsertBoundaryValuesEssential(dm, time, field, Nc, comps, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);CHKERRQ(ierr); 513 ierr = DMPlexLabelClearCells(dm,label);CHKERRQ(ierr); 514 break; 515 case DM_BC_ESSENTIAL_FIELD: 516 ierr = DMPlexLabelAddCells(dm,label);CHKERRQ(ierr); 517 ierr = DMPlexInsertBoundaryValuesEssentialField(dm, time, locX, field, Nc, comps, label, numids, ids, 518 (void (*)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 519 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 520 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) func, ctx, locX);CHKERRQ(ierr); 521 ierr = DMPlexLabelClearCells(dm,label);CHKERRQ(ierr); 522 break; 523 default: break; 524 } 525 } else if (id == PETSCFV_CLASSID) { 526 if (!faceGeomFVM) continue; 527 ierr = DMPlexInsertBoundaryValuesRiemann(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, field, Nc, comps, label, numids, ids, 528 (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);CHKERRQ(ierr); 529 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 530 } 531 PetscFunctionReturn(0); 532 } 533 534 /*@ 535 DMPlexInsertBoundaryValues - Puts coefficients which represent boundary values into the local solution vector 536 537 Input Parameters: 538 + dm - The DM 539 . insertEssential - Should I insert essential (e.g. Dirichlet) or inessential (e.g. Neumann) boundary conditions 540 . time - The time 541 . faceGeomFVM - Face geometry data for FV discretizations 542 . cellGeomFVM - Cell geometry data for FV discretizations 543 - gradFVM - Gradient reconstruction data for FV discretizations 544 545 Output Parameters: 546 . locX - Solution updated with boundary values 547 548 Level: developer 549 550 .seealso: DMProjectFunctionLabelLocal() 551 @*/ 552 PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM) 553 { 554 PetscErrorCode ierr; 555 556 PetscFunctionBegin; 557 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 558 PetscValidHeaderSpecific(locX, VEC_CLASSID, 2); 559 if (faceGeomFVM) {PetscValidHeaderSpecific(faceGeomFVM, VEC_CLASSID, 4);} 560 if (cellGeomFVM) {PetscValidHeaderSpecific(cellGeomFVM, VEC_CLASSID, 5);} 561 if (gradFVM) {PetscValidHeaderSpecific(gradFVM, VEC_CLASSID, 6);} 562 ierr = PetscTryMethod(dm,"DMPlexInsertBoundaryValues_C",(DM,PetscBool,Vec,PetscReal,Vec,Vec,Vec),(dm,insertEssential,locX,time,faceGeomFVM,cellGeomFVM,gradFVM));CHKERRQ(ierr); 563 PetscFunctionReturn(0); 564 } 565 566 PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 567 { 568 Vec localX; 569 PetscErrorCode ierr; 570 571 PetscFunctionBegin; 572 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 573 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, time, NULL, NULL, NULL);CHKERRQ(ierr); 574 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 575 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 576 ierr = DMPlexComputeL2DiffLocal(dm, time, funcs, ctxs, localX, diff);CHKERRQ(ierr); 577 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 578 PetscFunctionReturn(0); 579 } 580 581 /*@C 582 DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h. 583 584 Input Parameters: 585 + dm - The DM 586 . time - The time 587 . funcs - The functions to evaluate for each field component 588 . ctxs - Optional array of contexts to pass to each function, or NULL. 589 - localX - The coefficient vector u_h, a local vector 590 591 Output Parameter: 592 . diff - The diff ||u - u_h||_2 593 594 Level: developer 595 596 .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff() 597 @*/ 598 PetscErrorCode DMPlexComputeL2DiffLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec localX, PetscReal *diff) 599 { 600 const PetscInt debug = 0; 601 PetscSection section; 602 PetscQuadrature quad; 603 PetscScalar *funcVal, *interpolant; 604 PetscReal *coords, *detJ, *J; 605 PetscReal localDiff = 0.0; 606 const PetscReal *quadWeights; 607 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cellHeight, cStart, cEnd, cEndInterior, c, field, fieldOffset; 608 PetscErrorCode ierr; 609 610 PetscFunctionBegin; 611 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 612 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 613 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 614 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 615 for (field = 0; field < numFields; ++field) { 616 PetscObject obj; 617 PetscClassId id; 618 PetscInt Nc; 619 620 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 621 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 622 if (id == PETSCFE_CLASSID) { 623 PetscFE fe = (PetscFE) obj; 624 625 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 626 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 627 } else if (id == PETSCFV_CLASSID) { 628 PetscFV fv = (PetscFV) obj; 629 630 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 631 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 632 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 633 numComponents += Nc; 634 } 635 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);CHKERRQ(ierr); 636 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 637 ierr = PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);CHKERRQ(ierr); 638 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr); 639 ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr); 640 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 641 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 642 for (c = cStart; c < cEnd; ++c) { 643 PetscScalar *x = NULL; 644 PetscReal elemDiff = 0.0; 645 PetscInt qc = 0; 646 647 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);CHKERRQ(ierr); 648 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 649 650 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 651 PetscObject obj; 652 PetscClassId id; 653 void * const ctx = ctxs ? ctxs[field] : NULL; 654 PetscInt Nb, Nc, q, fc; 655 656 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 657 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 658 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 659 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 660 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 661 if (debug) { 662 char title[1024]; 663 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 664 ierr = DMPrintCellVector(c, title, Nb, &x[fieldOffset]);CHKERRQ(ierr); 665 } 666 for (q = 0; q < Nq; ++q) { 667 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); 668 ierr = (*funcs[field])(coordDim, time, &coords[coordDim * q], Nc, funcVal, ctx); 669 if (ierr) { 670 PetscErrorCode ierr2; 671 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 672 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 673 ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2); 674 CHKERRQ(ierr); 675 } 676 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 677 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 678 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 679 for (fc = 0; fc < Nc; ++fc) { 680 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 681 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);} 682 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]; 683 } 684 } 685 fieldOffset += Nb; 686 qc += Nc; 687 } 688 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 689 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 690 localDiff += elemDiff; 691 } 692 ierr = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr); 693 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 694 *diff = PetscSqrtReal(*diff); 695 PetscFunctionReturn(0); 696 } 697 698 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) 699 { 700 const PetscInt debug = 0; 701 PetscSection section; 702 PetscQuadrature quad; 703 Vec localX; 704 PetscScalar *funcVal, *interpolant; 705 const PetscReal *quadPoints, *quadWeights; 706 PetscReal *coords, *realSpaceDer, *J, *invJ, *detJ; 707 PetscReal localDiff = 0.0; 708 PetscInt dim, coordDim, qNc = 0, Nq = 0, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset; 709 PetscErrorCode ierr; 710 711 PetscFunctionBegin; 712 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 713 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 714 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 715 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 716 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 717 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 718 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 719 for (field = 0; field < numFields; ++field) { 720 PetscFE fe; 721 PetscInt Nc; 722 723 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 724 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 725 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 726 numComponents += Nc; 727 } 728 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 729 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 730 /* ierr = DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 731 ierr = PetscMalloc7(numComponents,&funcVal,coordDim*Nq,&coords,coordDim*Nq,&realSpaceDer,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ,numComponents,&interpolant,Nq,&detJ);CHKERRQ(ierr); 732 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 733 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 734 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 735 for (c = cStart; c < cEnd; ++c) { 736 PetscScalar *x = NULL; 737 PetscReal elemDiff = 0.0; 738 PetscInt qc = 0; 739 740 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, invJ, detJ);CHKERRQ(ierr); 741 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 742 743 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 744 PetscFE fe; 745 void * const ctx = ctxs ? ctxs[field] : NULL; 746 PetscReal *basisDer; 747 PetscInt Nb, Nc, q, fc; 748 749 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 750 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 751 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 752 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 753 if (debug) { 754 char title[1024]; 755 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 756 ierr = DMPrintCellVector(c, title, Nb, &x[fieldOffset]);CHKERRQ(ierr); 757 } 758 for (q = 0; q < Nq; ++q) { 759 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); 760 ierr = (*funcs[field])(coordDim, time, &coords[q*coordDim], n, numFields, funcVal, ctx); 761 if (ierr) { 762 PetscErrorCode ierr2; 763 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 764 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 765 ierr2 = PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);CHKERRQ(ierr2); 766 CHKERRQ(ierr); 767 } 768 ierr = PetscFEInterpolateGradient_Static(fe, &x[fieldOffset], coordDim, invJ, n, q, interpolant);CHKERRQ(ierr); 769 for (fc = 0; fc < Nc; ++fc) { 770 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 771 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);} 772 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]; 773 } 774 } 775 fieldOffset += Nb; 776 qc += Nc; 777 } 778 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 779 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 780 localDiff += elemDiff; 781 } 782 ierr = PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);CHKERRQ(ierr); 783 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 784 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 785 *diff = PetscSqrtReal(*diff); 786 PetscFunctionReturn(0); 787 } 788 789 PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 790 { 791 const PetscInt debug = 0; 792 PetscSection section; 793 PetscQuadrature quad; 794 Vec localX; 795 PetscScalar *funcVal, *interpolant; 796 PetscReal *coords, *detJ, *J; 797 PetscReal *localDiff; 798 const PetscReal *quadPoints, *quadWeights; 799 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset; 800 PetscErrorCode ierr; 801 802 PetscFunctionBegin; 803 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 804 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 805 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 806 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 807 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 808 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 809 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 810 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 811 for (field = 0; field < numFields; ++field) { 812 PetscObject obj; 813 PetscClassId id; 814 PetscInt Nc; 815 816 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 817 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 818 if (id == PETSCFE_CLASSID) { 819 PetscFE fe = (PetscFE) obj; 820 821 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 822 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 823 } else if (id == PETSCFV_CLASSID) { 824 PetscFV fv = (PetscFV) obj; 825 826 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 827 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 828 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 829 numComponents += Nc; 830 } 831 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 832 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 833 ierr = PetscCalloc6(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);CHKERRQ(ierr); 834 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 835 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 836 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 837 for (c = cStart; c < cEnd; ++c) { 838 PetscScalar *x = NULL; 839 PetscInt qc = 0; 840 841 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);CHKERRQ(ierr); 842 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 843 844 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 845 PetscObject obj; 846 PetscClassId id; 847 void * const ctx = ctxs ? ctxs[field] : NULL; 848 PetscInt Nb, Nc, q, fc; 849 850 PetscReal elemDiff = 0.0; 851 852 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 853 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 854 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 855 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 856 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 857 if (debug) { 858 char title[1024]; 859 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 860 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 861 } 862 for (q = 0; q < Nq; ++q) { 863 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); 864 ierr = (*funcs[field])(coordDim, time, &coords[coordDim*q], numFields, funcVal, ctx); 865 if (ierr) { 866 PetscErrorCode ierr2; 867 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 868 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 869 ierr2 = PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2); 870 CHKERRQ(ierr); 871 } 872 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 873 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 874 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 875 for (fc = 0; fc < Nc; ++fc) { 876 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 877 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);} 878 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]; 879 } 880 } 881 fieldOffset += Nb; 882 qc += Nc; 883 localDiff[field] += elemDiff; 884 } 885 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 886 } 887 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 888 ierr = MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 889 for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]); 890 ierr = PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr); 891 PetscFunctionReturn(0); 892 } 893 894 /*@C 895 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. 896 897 Input Parameters: 898 + dm - The DM 899 . time - The time 900 . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation 901 . ctxs - Optional array of contexts to pass to each function, or NULL. 902 - X - The coefficient vector u_h 903 904 Output Parameter: 905 . D - A Vec which holds the difference ||u - u_h||_2 for each cell 906 907 Level: developer 908 909 .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff() 910 @*/ 911 PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D) 912 { 913 PetscSection section; 914 PetscQuadrature quad; 915 Vec localX; 916 PetscScalar *funcVal, *interpolant; 917 PetscReal *coords, *detJ, *J; 918 const PetscReal *quadPoints, *quadWeights; 919 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset; 920 PetscErrorCode ierr; 921 922 PetscFunctionBegin; 923 ierr = VecSet(D, 0.0);CHKERRQ(ierr); 924 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 925 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 926 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 927 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 928 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 929 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 930 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 931 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 932 for (field = 0; field < numFields; ++field) { 933 PetscObject obj; 934 PetscClassId id; 935 PetscInt Nc; 936 937 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 938 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 939 if (id == PETSCFE_CLASSID) { 940 PetscFE fe = (PetscFE) obj; 941 942 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 943 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 944 } else if (id == PETSCFV_CLASSID) { 945 PetscFV fv = (PetscFV) obj; 946 947 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 948 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 949 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 950 numComponents += Nc; 951 } 952 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 953 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 954 ierr = PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);CHKERRQ(ierr); 955 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 956 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 957 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 958 for (c = cStart; c < cEnd; ++c) { 959 PetscScalar *x = NULL; 960 PetscScalar elemDiff = 0.0; 961 PetscInt qc = 0; 962 963 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);CHKERRQ(ierr); 964 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 965 966 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 967 PetscObject obj; 968 PetscClassId id; 969 void * const ctx = ctxs ? ctxs[field] : NULL; 970 PetscInt Nb, Nc, q, fc; 971 972 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 973 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 974 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 975 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 976 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 977 if (funcs[field]) { 978 for (q = 0; q < Nq; ++q) { 979 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); 980 ierr = (*funcs[field])(dim, time, &coords[q*coordDim], Nc, funcVal, ctx); 981 if (ierr) { 982 PetscErrorCode ierr2; 983 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 984 ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2); 985 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 986 CHKERRQ(ierr); 987 } 988 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 989 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 990 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 991 for (fc = 0; fc < Nc; ++fc) { 992 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 993 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]; 994 } 995 } 996 } 997 fieldOffset += Nb; 998 qc += Nc; 999 } 1000 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1001 ierr = VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);CHKERRQ(ierr); 1002 } 1003 ierr = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr); 1004 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1005 ierr = VecSqrtAbs(D);CHKERRQ(ierr); 1006 PetscFunctionReturn(0); 1007 } 1008 1009 /*@C 1010 DMPlexComputeGradientClementInterpolant - This function computes the L2 projection of the cellwise gradient of a function u onto P1, and stores it in a Vec. 1011 1012 Input Parameters: 1013 + dm - The DM 1014 - LocX - The coefficient vector u_h 1015 1016 Output Parameter: 1017 . locC - A Vec which holds the Clement interpolant of the gradient 1018 1019 Notes: Add citation to (Clement, 1975) and definition of the interpolant 1020 \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 1021 1022 Level: developer 1023 1024 .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff() 1025 @*/ 1026 PetscErrorCode DMPlexComputeGradientClementInterpolant(DM dm, Vec locX, Vec locC) 1027 { 1028 DM_Plex *mesh = (DM_Plex *) dm->data; 1029 PetscInt debug = mesh->printFEM; 1030 DM dmC; 1031 PetscSection section; 1032 PetscQuadrature quad; 1033 PetscScalar *interpolant, *gradsum; 1034 PetscReal *coords, *detJ, *J, *invJ; 1035 const PetscReal *quadPoints, *quadWeights; 1036 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, vStart, vEnd, v, field, fieldOffset; 1037 PetscErrorCode ierr; 1038 1039 PetscFunctionBegin; 1040 ierr = VecGetDM(locC, &dmC);CHKERRQ(ierr); 1041 ierr = VecSet(locC, 0.0);CHKERRQ(ierr); 1042 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1043 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 1044 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1045 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1046 for (field = 0; field < numFields; ++field) { 1047 PetscObject obj; 1048 PetscClassId id; 1049 PetscInt Nc; 1050 1051 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1052 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1053 if (id == PETSCFE_CLASSID) { 1054 PetscFE fe = (PetscFE) obj; 1055 1056 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1057 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1058 } else if (id == PETSCFV_CLASSID) { 1059 PetscFV fv = (PetscFV) obj; 1060 1061 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1062 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1063 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1064 numComponents += Nc; 1065 } 1066 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 1067 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 1068 ierr = PetscMalloc6(coordDim*numComponents*2,&gradsum,coordDim*numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ);CHKERRQ(ierr); 1069 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1070 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1071 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1072 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1073 for (v = vStart; v < vEnd; ++v) { 1074 PetscScalar volsum = 0.0; 1075 PetscInt *star = NULL; 1076 PetscInt starSize, st, d, fc; 1077 1078 ierr = PetscMemzero(gradsum, coordDim*numComponents * sizeof(PetscScalar));CHKERRQ(ierr); 1079 ierr = DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1080 for (st = 0; st < starSize*2; st += 2) { 1081 const PetscInt cell = star[st]; 1082 PetscScalar *grad = &gradsum[coordDim*numComponents]; 1083 PetscScalar *x = NULL; 1084 PetscReal vol = 0.0; 1085 1086 if ((cell < cStart) || (cell >= cEnd)) continue; 1087 ierr = DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, J, invJ, detJ);CHKERRQ(ierr); 1088 ierr = DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr); 1089 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1090 PetscObject obj; 1091 PetscClassId id; 1092 PetscInt Nb, Nc, q, qc = 0; 1093 1094 ierr = PetscMemzero(grad, coordDim*numComponents * sizeof(PetscScalar));CHKERRQ(ierr); 1095 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1096 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1097 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1098 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1099 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1100 for (q = 0; q < Nq; ++q) { 1101 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); 1102 if (ierr) { 1103 PetscErrorCode ierr2; 1104 ierr2 = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr2); 1105 ierr2 = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr2); 1106 ierr2 = PetscFree4(interpolant,coords,detJ,J);CHKERRQ(ierr2); 1107 CHKERRQ(ierr); 1108 } 1109 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolateGradient_Static((PetscFE) obj, &x[fieldOffset], coordDim, invJ, NULL, q, interpolant);CHKERRQ(ierr);} 1110 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1111 for (fc = 0; fc < Nc; ++fc) { 1112 const PetscReal wt = quadWeights[q*qNc+qc+fc]; 1113 1114 for (d = 0; d < coordDim; ++d) grad[fc*coordDim+d] += interpolant[fc*dim+d]*wt*detJ[q]; 1115 } 1116 vol += quadWeights[q*qNc]*detJ[q]; 1117 } 1118 fieldOffset += Nb; 1119 qc += Nc; 1120 } 1121 ierr = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr); 1122 for (fc = 0; fc < numComponents; ++fc) { 1123 for (d = 0; d < coordDim; ++d) { 1124 gradsum[fc*coordDim+d] += grad[fc*coordDim+d]; 1125 } 1126 } 1127 volsum += vol; 1128 if (debug) { 1129 ierr = PetscPrintf(PETSC_COMM_SELF, "Cell %D gradient: [", cell);CHKERRQ(ierr); 1130 for (fc = 0; fc < numComponents; ++fc) { 1131 for (d = 0; d < coordDim; ++d) { 1132 if (fc || d > 0) {ierr = PetscPrintf(PETSC_COMM_SELF, ", ");CHKERRQ(ierr);} 1133 ierr = PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc*coordDim+d]));CHKERRQ(ierr); 1134 } 1135 } 1136 ierr = PetscPrintf(PETSC_COMM_SELF, "]\n");CHKERRQ(ierr); 1137 } 1138 } 1139 for (fc = 0; fc < numComponents; ++fc) { 1140 for (d = 0; d < coordDim; ++d) gradsum[fc*coordDim+d] /= volsum; 1141 } 1142 ierr = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1143 ierr = DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES);CHKERRQ(ierr); 1144 } 1145 ierr = PetscFree6(gradsum,interpolant,coords,detJ,J,invJ);CHKERRQ(ierr); 1146 PetscFunctionReturn(0); 1147 } 1148 1149 static PetscErrorCode DMPlexComputeIntegral_Internal(DM dm, Vec X, PetscInt cStart, PetscInt cEnd, PetscScalar *cintegral, void *user) 1150 { 1151 DM dmAux = NULL; 1152 PetscDS prob, probAux; 1153 PetscSection section, sectionAux; 1154 Vec locX, locA; 1155 PetscInt dim, numCells = cEnd - cStart, cell, c, f; 1156 PetscBool useFEM = PETSC_FALSE, useFVM = PETSC_FALSE; 1157 /* DS */ 1158 PetscInt Nf, totDim, *uOff, *uOff_x, numConstants; 1159 PetscInt NfAux, totDimAux, *aOff; 1160 PetscScalar *u, *a; 1161 const PetscScalar *constants; 1162 /* Geometry */ 1163 PetscFECellGeom *cgeomFEM; 1164 DM dmGrad; 1165 Vec cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL; 1166 PetscFVCellGeom *cgeomFVM; 1167 const PetscScalar *lgrad; 1168 PetscErrorCode ierr; 1169 1170 PetscFunctionBegin; 1171 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1172 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1173 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1174 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1175 /* Determine which discretizations we have */ 1176 for (f = 0; f < Nf; ++f) { 1177 PetscObject obj; 1178 PetscClassId id; 1179 1180 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1181 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1182 if (id == PETSCFE_CLASSID) useFEM = PETSC_TRUE; 1183 if (id == PETSCFV_CLASSID) useFVM = PETSC_TRUE; 1184 } 1185 /* Get local solution with boundary values */ 1186 ierr = DMGetLocalVector(dm, &locX);CHKERRQ(ierr); 1187 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 1188 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);CHKERRQ(ierr); 1189 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);CHKERRQ(ierr); 1190 /* Read DS information */ 1191 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1192 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 1193 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 1194 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 1195 /* Read Auxiliary DS information */ 1196 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1197 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);CHKERRQ(ierr); 1198 if (dmAux) { 1199 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1200 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 1201 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1202 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1203 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 1204 } 1205 /* Allocate data arrays */ 1206 ierr = PetscCalloc1(numCells*totDim, &u);CHKERRQ(ierr); 1207 if (useFEM) {ierr = PetscMalloc1(numCells, &cgeomFEM);CHKERRQ(ierr);} 1208 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1209 /* Read out geometry */ 1210 if (useFEM) { 1211 for (cell = cStart; cell < cEnd; ++cell) { 1212 const PetscInt c = cell - cStart; 1213 1214 ierr = DMPlexComputeCellGeometryFEM(dm, cell, NULL, cgeomFEM[c].v0, cgeomFEM[c].J, cgeomFEM[c].invJ, &cgeomFEM[c].detJ);CHKERRQ(ierr); 1215 if (cgeomFEM[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeomFEM[c].detJ, c); 1216 } 1217 } 1218 if (useFVM) { 1219 PetscFV fv = NULL; 1220 Vec grad; 1221 PetscInt fStart, fEnd; 1222 PetscBool compGrad; 1223 1224 for (f = 0; f < Nf; ++f) { 1225 PetscObject obj; 1226 PetscClassId id; 1227 1228 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1229 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1230 if (id == PETSCFV_CLASSID) {fv = (PetscFV) obj; break;} 1231 } 1232 ierr = PetscFVGetComputeGradients(fv, &compGrad);CHKERRQ(ierr); 1233 ierr = PetscFVSetComputeGradients(fv, PETSC_TRUE);CHKERRQ(ierr); 1234 ierr = DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);CHKERRQ(ierr); 1235 ierr = DMPlexComputeGradientFVM(dm, fv, faceGeometryFVM, cellGeometryFVM, &dmGrad);CHKERRQ(ierr); 1236 ierr = PetscFVSetComputeGradients(fv, compGrad);CHKERRQ(ierr); 1237 ierr = VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1238 /* Reconstruct and limit cell gradients */ 1239 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1240 ierr = DMGetGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1241 ierr = DMPlexReconstructGradients_Internal(dm, fv, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);CHKERRQ(ierr); 1242 /* Communicate gradient values */ 1243 ierr = DMGetLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1244 ierr = DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1245 ierr = DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1246 ierr = DMRestoreGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1247 /* Handle non-essential (e.g. outflow) boundary values */ 1248 ierr = DMPlexInsertBoundaryValues(dm, PETSC_FALSE, locX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);CHKERRQ(ierr); 1249 ierr = VecGetArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1250 } 1251 /* Read out data from inputs */ 1252 for (c = cStart; c < cEnd; ++c) { 1253 PetscScalar *x = NULL; 1254 PetscInt i; 1255 1256 ierr = DMPlexVecGetClosure(dm, section, locX, c, NULL, &x);CHKERRQ(ierr); 1257 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 1258 ierr = DMPlexVecRestoreClosure(dm, section, locX, c, NULL, &x);CHKERRQ(ierr); 1259 if (dmAux) { 1260 ierr = DMPlexVecGetClosure(dmAux, sectionAux, locA, c, NULL, &x);CHKERRQ(ierr); 1261 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 1262 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, locA, c, NULL, &x);CHKERRQ(ierr); 1263 } 1264 } 1265 /* Do integration for each field */ 1266 for (f = 0; f < Nf; ++f) { 1267 PetscObject obj; 1268 PetscClassId id; 1269 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 1270 1271 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1272 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1273 if (id == PETSCFE_CLASSID) { 1274 PetscFE fe = (PetscFE) obj; 1275 PetscQuadrature q; 1276 PetscInt Nq, Nb; 1277 1278 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1279 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1280 ierr = PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 1281 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1282 blockSize = Nb*Nq; 1283 batchSize = numBlocks * blockSize; 1284 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1285 numChunks = numCells / (numBatches*batchSize); 1286 Ne = numChunks*numBatches*batchSize; 1287 Nr = numCells % (numBatches*batchSize); 1288 offset = numCells - Nr; 1289 ierr = PetscFEIntegrate(fe, prob, f, Ne, cgeomFEM, u, probAux, a, cintegral);CHKERRQ(ierr); 1290 ierr = PetscFEIntegrate(fe, prob, f, Nr, &cgeomFEM[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], &cintegral[offset*Nf]);CHKERRQ(ierr); 1291 } else if (id == PETSCFV_CLASSID) { 1292 PetscInt foff; 1293 PetscPointFunc obj_func; 1294 PetscScalar lint; 1295 1296 ierr = PetscDSGetObjective(prob, f, &obj_func);CHKERRQ(ierr); 1297 ierr = PetscDSGetFieldOffset(prob, f, &foff);CHKERRQ(ierr); 1298 if (obj_func) { 1299 for (c = 0; c < numCells; ++c) { 1300 PetscScalar *u_x; 1301 1302 ierr = DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);CHKERRQ(ierr); 1303 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); 1304 cintegral[c*Nf+f] += PetscRealPart(lint)*cgeomFVM[c].volume; 1305 } 1306 } 1307 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 1308 } 1309 /* Cleanup data arrays */ 1310 if (useFVM) { 1311 ierr = VecRestoreArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1312 ierr = VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1313 ierr = DMRestoreLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1314 ierr = VecDestroy(&faceGeometryFVM);CHKERRQ(ierr); 1315 ierr = VecDestroy(&cellGeometryFVM);CHKERRQ(ierr); 1316 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr); 1317 } 1318 if (useFEM) {ierr = PetscFree(cgeomFEM);CHKERRQ(ierr);} 1319 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1320 ierr = PetscFree(u);CHKERRQ(ierr); 1321 /* Cleanup */ 1322 ierr = DMRestoreLocalVector(dm, &locX);CHKERRQ(ierr); 1323 PetscFunctionReturn(0); 1324 } 1325 1326 /*@ 1327 DMPlexComputeIntegralFEM - Form the integral over the domain from the global input X using pointwise functions specified by the user 1328 1329 Input Parameters: 1330 + dm - The mesh 1331 . X - Global input vector 1332 - user - The user context 1333 1334 Output Parameter: 1335 . integral - Integral for each field 1336 1337 Level: developer 1338 1339 .seealso: DMPlexComputeResidualFEM() 1340 @*/ 1341 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user) 1342 { 1343 DM_Plex *mesh = (DM_Plex *) dm->data; 1344 PetscScalar *cintegral; 1345 PetscReal *lintegral; 1346 PetscInt Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell; 1347 PetscErrorCode ierr; 1348 1349 PetscFunctionBegin; 1350 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1351 PetscValidHeaderSpecific(X, VEC_CLASSID, 2); 1352 PetscValidPointer(integral, 3); 1353 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1354 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 1355 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr); 1356 ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr); 1357 ierr = DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);CHKERRQ(ierr); 1358 cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight]; 1359 /* TODO Introduce a loop over large chunks (right now this is a single chunk) */ 1360 ierr = PetscCalloc2(Nf, &lintegral, (cEnd-cStart)*Nf, &cintegral);CHKERRQ(ierr); 1361 ierr = DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);CHKERRQ(ierr); 1362 /* Sum up values */ 1363 for (cell = cStart; cell < cEnd; ++cell) { 1364 const PetscInt c = cell - cStart; 1365 1366 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);CHKERRQ(ierr);} 1367 for (f = 0; f < Nf; ++f) lintegral[f] += PetscRealPart(cintegral[c*Nf+f]); 1368 } 1369 ierr = MPIU_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr); 1370 if (mesh->printFEM) { 1371 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Integral:");CHKERRQ(ierr); 1372 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", integral[f]);CHKERRQ(ierr);} 1373 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 1374 } 1375 ierr = PetscFree2(lintegral, cintegral);CHKERRQ(ierr); 1376 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1377 PetscFunctionReturn(0); 1378 } 1379 1380 /*@ 1381 DMPlexComputeCellwiseIntegralFEM - Form the vector of cellwise integrals F from the global input X using pointwise functions specified by the user 1382 1383 Input Parameters: 1384 + dm - The mesh 1385 . X - Global input vector 1386 - user - The user context 1387 1388 Output Parameter: 1389 . integral - Cellwise integrals for each field 1390 1391 Level: developer 1392 1393 .seealso: DMPlexComputeResidualFEM() 1394 @*/ 1395 PetscErrorCode DMPlexComputeCellwiseIntegralFEM(DM dm, Vec X, Vec F, void *user) 1396 { 1397 DM_Plex *mesh = (DM_Plex *) dm->data; 1398 DM dmF; 1399 PetscSection sectionF; 1400 PetscScalar *cintegral, *af; 1401 PetscInt Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell; 1402 PetscErrorCode ierr; 1403 1404 PetscFunctionBegin; 1405 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1406 PetscValidHeaderSpecific(X, VEC_CLASSID, 2); 1407 PetscValidHeaderSpecific(F, VEC_CLASSID, 3); 1408 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1409 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 1410 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr); 1411 ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr); 1412 ierr = DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);CHKERRQ(ierr); 1413 cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight]; 1414 /* TODO Introduce a loop over large chunks (right now this is a single chunk) */ 1415 ierr = PetscCalloc1((cEnd-cStart)*Nf, &cintegral);CHKERRQ(ierr); 1416 ierr = DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);CHKERRQ(ierr); 1417 /* Put values in F*/ 1418 ierr = VecGetDM(F, &dmF);CHKERRQ(ierr); 1419 ierr = DMGetDefaultSection(dmF, §ionF);CHKERRQ(ierr); 1420 ierr = VecGetArray(F, &af);CHKERRQ(ierr); 1421 for (cell = cStart; cell < cEnd; ++cell) { 1422 const PetscInt c = cell - cStart; 1423 PetscInt dof, off; 1424 1425 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);CHKERRQ(ierr);} 1426 ierr = PetscSectionGetDof(sectionF, cell, &dof);CHKERRQ(ierr); 1427 ierr = PetscSectionGetOffset(sectionF, cell, &off);CHKERRQ(ierr); 1428 if (dof != Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cell dofs %D != %D", dof, Nf); 1429 for (f = 0; f < Nf; ++f) af[off+f] = cintegral[c*Nf+f]; 1430 } 1431 ierr = VecRestoreArray(F, &af);CHKERRQ(ierr); 1432 ierr = PetscFree(cintegral);CHKERRQ(ierr); 1433 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1434 PetscFunctionReturn(0); 1435 } 1436 1437 /*@ 1438 DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1439 1440 Input Parameters: 1441 + dmf - The fine mesh 1442 . dmc - The coarse mesh 1443 - user - The user context 1444 1445 Output Parameter: 1446 . In - The interpolation matrix 1447 1448 Level: developer 1449 1450 .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1451 @*/ 1452 PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user) 1453 { 1454 DM_Plex *mesh = (DM_Plex *) dmc->data; 1455 const char *name = "Interpolator"; 1456 PetscDS prob; 1457 PetscFE *feRef; 1458 PetscFV *fvRef; 1459 PetscSection fsection, fglobalSection; 1460 PetscSection csection, cglobalSection; 1461 PetscScalar *elemMat; 1462 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c; 1463 PetscInt cTotDim, rTotDim = 0; 1464 PetscErrorCode ierr; 1465 1466 PetscFunctionBegin; 1467 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1468 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1469 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1470 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1471 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1472 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1473 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1474 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1475 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1476 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1477 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr); 1478 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1479 for (f = 0; f < Nf; ++f) { 1480 PetscObject obj; 1481 PetscClassId id; 1482 PetscInt rNb = 0, Nc = 0; 1483 1484 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1485 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1486 if (id == PETSCFE_CLASSID) { 1487 PetscFE fe = (PetscFE) obj; 1488 1489 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1490 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1491 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1492 } else if (id == PETSCFV_CLASSID) { 1493 PetscFV fv = (PetscFV) obj; 1494 PetscDualSpace Q; 1495 1496 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1497 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1498 ierr = PetscDualSpaceGetDimension(Q, &rNb);CHKERRQ(ierr); 1499 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1500 } 1501 rTotDim += rNb; 1502 } 1503 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1504 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 1505 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 1506 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1507 PetscDualSpace Qref; 1508 PetscQuadrature f; 1509 const PetscReal *qpoints, *qweights; 1510 PetscReal *points; 1511 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1512 1513 /* Compose points from all dual basis functionals */ 1514 if (feRef[fieldI]) { 1515 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1516 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 1517 } else { 1518 ierr = PetscFVGetDualSpace(fvRef[fieldI], &Qref);CHKERRQ(ierr); 1519 ierr = PetscFVGetNumComponents(fvRef[fieldI], &Nc);CHKERRQ(ierr); 1520 } 1521 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1522 for (i = 0; i < fpdim; ++i) { 1523 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1524 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1525 npoints += Np; 1526 } 1527 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1528 for (i = 0, k = 0; i < fpdim; ++i) { 1529 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1530 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1531 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1532 } 1533 1534 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1535 PetscObject obj; 1536 PetscClassId id; 1537 PetscReal *B; 1538 PetscInt NcJ = 0, cpdim = 0, j, qNc; 1539 1540 ierr = PetscDSGetDiscretization(prob, fieldJ, &obj);CHKERRQ(ierr); 1541 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1542 if (id == PETSCFE_CLASSID) { 1543 PetscFE fe = (PetscFE) obj; 1544 1545 /* Evaluate basis at points */ 1546 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 1547 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1548 /* For now, fields only interpolate themselves */ 1549 if (fieldI == fieldJ) { 1550 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); 1551 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1552 for (i = 0, k = 0; i < fpdim; ++i) { 1553 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1554 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);CHKERRQ(ierr); 1555 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); 1556 for (p = 0; p < Np; ++p, ++k) { 1557 for (j = 0; j < cpdim; ++j) { 1558 /* NOTE: This is not quite right, unless fpdim == number of fine grid functional quad points */ 1559 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p*qNc+c]; 1560 } 1561 } 1562 } 1563 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1564 } 1565 } else if (id == PETSCFV_CLASSID) { 1566 PetscFV fv = (PetscFV) obj; 1567 1568 /* Evaluate constant function at points */ 1569 ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr); 1570 cpdim = 1; 1571 /* For now, fields only interpolate themselves */ 1572 if (fieldI == fieldJ) { 1573 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); 1574 for (i = 0, k = 0; i < fpdim; ++i) { 1575 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1576 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);CHKERRQ(ierr); 1577 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); 1578 for (p = 0; p < Np; ++p, ++k) { 1579 for (j = 0; j < cpdim; ++j) { 1580 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p*qNc+c]; 1581 } 1582 } 1583 } 1584 } 1585 } 1586 offsetJ += cpdim*NcJ; 1587 } 1588 offsetI += fpdim*Nc; 1589 ierr = PetscFree(points);CHKERRQ(ierr); 1590 } 1591 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1592 /* Preallocate matrix */ 1593 { 1594 Mat preallocator; 1595 PetscScalar *vals; 1596 PetscInt *cellCIndices, *cellFIndices; 1597 PetscInt locRows, locCols, cell; 1598 1599 ierr = MatGetLocalSize(In, &locRows, &locCols);CHKERRQ(ierr); 1600 ierr = MatCreate(PetscObjectComm((PetscObject) In), &preallocator);CHKERRQ(ierr); 1601 ierr = MatSetType(preallocator, MATPREALLOCATOR);CHKERRQ(ierr); 1602 ierr = MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1603 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 1604 ierr = PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1605 for (cell = cStart; cell < cEnd; ++cell) { 1606 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1607 ierr = MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);CHKERRQ(ierr); 1608 } 1609 ierr = PetscFree3(vals,cellCIndices,cellFIndices);CHKERRQ(ierr); 1610 ierr = MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1611 ierr = MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1612 ierr = MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);CHKERRQ(ierr); 1613 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 1614 } 1615 /* Fill matrix */ 1616 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1617 for (c = cStart; c < cEnd; ++c) { 1618 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1619 } 1620 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1621 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1622 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1623 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1624 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1625 if (mesh->printFEM) { 1626 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1627 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1628 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1629 } 1630 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1631 PetscFunctionReturn(0); 1632 } 1633 1634 PetscErrorCode DMPlexComputeMassMatrixNested(DM dmc, DM dmf, Mat mass, void *user) 1635 { 1636 SETERRQ(PetscObjectComm((PetscObject) dmc), PETSC_ERR_SUP, "Laziness"); 1637 } 1638 1639 /*@ 1640 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 1641 1642 Input Parameters: 1643 + dmf - The fine mesh 1644 . dmc - The coarse mesh 1645 - user - The user context 1646 1647 Output Parameter: 1648 . In - The interpolation matrix 1649 1650 Level: developer 1651 1652 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1653 @*/ 1654 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 1655 { 1656 DM_Plex *mesh = (DM_Plex *) dmf->data; 1657 const char *name = "Interpolator"; 1658 PetscDS prob; 1659 PetscSection fsection, csection, globalFSection, globalCSection; 1660 PetscHashJK ht; 1661 PetscLayout rLayout; 1662 PetscInt *dnz, *onz; 1663 PetscInt locRows, rStart, rEnd; 1664 PetscReal *x, *v0, *J, *invJ, detJ; 1665 PetscReal *v0c, *Jc, *invJc, detJc; 1666 PetscScalar *elemMat; 1667 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1668 PetscErrorCode ierr; 1669 1670 PetscFunctionBegin; 1671 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1672 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1673 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1674 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1675 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1676 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1677 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1678 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1679 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1680 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1681 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1682 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1683 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1684 ierr = PetscMalloc1(totDim, &elemMat);CHKERRQ(ierr); 1685 1686 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1687 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1688 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1689 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1690 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1691 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1692 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1693 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1694 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1695 for (field = 0; field < Nf; ++field) { 1696 PetscObject obj; 1697 PetscClassId id; 1698 PetscDualSpace Q = NULL; 1699 PetscQuadrature f; 1700 const PetscReal *qpoints; 1701 PetscInt Nc, Np, fpdim, i, d; 1702 1703 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1704 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1705 if (id == PETSCFE_CLASSID) { 1706 PetscFE fe = (PetscFE) obj; 1707 1708 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1709 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1710 } else if (id == PETSCFV_CLASSID) { 1711 PetscFV fv = (PetscFV) obj; 1712 1713 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1714 Nc = 1; 1715 } 1716 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1717 /* For each fine grid cell */ 1718 for (cell = cStart; cell < cEnd; ++cell) { 1719 PetscInt *findices, *cindices; 1720 PetscInt numFIndices, numCIndices; 1721 1722 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1723 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1724 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim); 1725 for (i = 0; i < fpdim; ++i) { 1726 Vec pointVec; 1727 PetscScalar *pV; 1728 PetscSF coarseCellSF = NULL; 1729 const PetscSFNode *coarseCells; 1730 PetscInt numCoarseCells, q, c; 1731 1732 /* Get points from the dual basis functional quadrature */ 1733 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1734 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1735 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1736 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1737 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1738 for (q = 0; q < Np; ++q) { 1739 /* Transform point to real space */ 1740 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1741 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1742 } 1743 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1744 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1745 /* OPT: Pack all quad points from fine cell */ 1746 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1747 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 1748 /* Update preallocation info */ 1749 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1750 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1751 { 1752 PetscHashJKKey key; 1753 PetscHashJKIter missing, iter; 1754 1755 key.j = findices[i]; 1756 if (key.j >= 0) { 1757 /* Get indices for coarse elements */ 1758 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1759 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1760 for (c = 0; c < numCIndices; ++c) { 1761 key.k = cindices[c]; 1762 if (key.k < 0) continue; 1763 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1764 if (missing) { 1765 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1766 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1767 else ++onz[key.j-rStart]; 1768 } 1769 } 1770 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1771 } 1772 } 1773 } 1774 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1775 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1776 } 1777 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1778 } 1779 } 1780 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1781 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1782 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1783 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1784 for (field = 0; field < Nf; ++field) { 1785 PetscObject obj; 1786 PetscClassId id; 1787 PetscDualSpace Q = NULL; 1788 PetscQuadrature f; 1789 const PetscReal *qpoints, *qweights; 1790 PetscInt Nc, qNc, Np, fpdim, i, d; 1791 1792 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1793 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1794 if (id == PETSCFE_CLASSID) { 1795 PetscFE fe = (PetscFE) obj; 1796 1797 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1798 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1799 } else if (id == PETSCFV_CLASSID) { 1800 PetscFV fv = (PetscFV) obj; 1801 1802 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1803 Nc = 1; 1804 } 1805 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1806 /* For each fine grid cell */ 1807 for (cell = cStart; cell < cEnd; ++cell) { 1808 PetscInt *findices, *cindices; 1809 PetscInt numFIndices, numCIndices; 1810 1811 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1812 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1813 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim); 1814 for (i = 0; i < fpdim; ++i) { 1815 Vec pointVec; 1816 PetscScalar *pV; 1817 PetscSF coarseCellSF = NULL; 1818 const PetscSFNode *coarseCells; 1819 PetscInt numCoarseCells, cpdim, q, c, j; 1820 1821 /* Get points from the dual basis functional quadrature */ 1822 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1823 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);CHKERRQ(ierr); 1824 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); 1825 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1826 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1827 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1828 for (q = 0; q < Np; ++q) { 1829 /* Transform point to real space */ 1830 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1831 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1832 } 1833 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1834 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1835 /* OPT: Read this out from preallocation information */ 1836 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1837 /* Update preallocation info */ 1838 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1839 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1840 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1841 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1842 PetscReal pVReal[3]; 1843 1844 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1845 /* Transform points from real space to coarse reference space */ 1846 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 1847 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 1848 CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x); 1849 1850 if (id == PETSCFE_CLASSID) { 1851 PetscFE fe = (PetscFE) obj; 1852 PetscReal *B; 1853 1854 /* Evaluate coarse basis on contained point */ 1855 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1856 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 1857 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 1858 /* Get elemMat entries by multiplying by weight */ 1859 for (j = 0; j < cpdim; ++j) { 1860 for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*qweights[ccell*qNc + c]; 1861 } 1862 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1863 } else { 1864 cpdim = 1; 1865 for (j = 0; j < cpdim; ++j) { 1866 for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c]; 1867 } 1868 } 1869 /* Update interpolator */ 1870 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 1871 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 1872 ierr = MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1873 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1874 } 1875 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1876 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1877 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1878 } 1879 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1880 } 1881 } 1882 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 1883 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 1884 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1885 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1886 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1887 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1888 PetscFunctionReturn(0); 1889 } 1890 1891 /*@ 1892 DMPlexComputeMassMatrixGeneral - Form the local portion of the mass matrix M from the coarse DM to a non-nested fine DM. 1893 1894 Input Parameters: 1895 + dmf - The fine mesh 1896 . dmc - The coarse mesh 1897 - user - The user context 1898 1899 Output Parameter: 1900 . mass - The mass matrix 1901 1902 Level: developer 1903 1904 .seealso: DMPlexComputeMassMatrixNested(), DMPlexComputeInterpolatorNested(), DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1905 @*/ 1906 PetscErrorCode DMPlexComputeMassMatrixGeneral(DM dmc, DM dmf, Mat mass, void *user) 1907 { 1908 DM_Plex *mesh = (DM_Plex *) dmf->data; 1909 const char *name = "Mass Matrix"; 1910 PetscDS prob; 1911 PetscSection fsection, csection, globalFSection, globalCSection; 1912 PetscHashJK ht; 1913 PetscLayout rLayout; 1914 PetscInt *dnz, *onz; 1915 PetscInt locRows, rStart, rEnd; 1916 PetscReal *x, *v0, *J, *invJ, detJ; 1917 PetscReal *v0c, *Jc, *invJc, detJc; 1918 PetscScalar *elemMat; 1919 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1920 PetscErrorCode ierr; 1921 1922 PetscFunctionBegin; 1923 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1924 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1925 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1926 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1927 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1928 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1929 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1930 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1931 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1932 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1933 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1934 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1935 ierr = PetscMalloc1(totDim, &elemMat);CHKERRQ(ierr); 1936 1937 ierr = MatGetLocalSize(mass, &locRows, NULL);CHKERRQ(ierr); 1938 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) mass), &rLayout);CHKERRQ(ierr); 1939 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1940 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1941 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1942 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1943 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1944 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1945 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1946 for (field = 0; field < Nf; ++field) { 1947 PetscObject obj; 1948 PetscClassId id; 1949 PetscQuadrature quad; 1950 const PetscReal *qpoints; 1951 PetscInt Nq, Nc, i, d; 1952 1953 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1954 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1955 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetQuadrature((PetscFE) obj, &quad);CHKERRQ(ierr);} 1956 else {ierr = PetscFVGetQuadrature((PetscFV) obj, &quad);CHKERRQ(ierr);} 1957 ierr = PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, NULL);CHKERRQ(ierr); 1958 /* For each fine grid cell */ 1959 for (cell = cStart; cell < cEnd; ++cell) { 1960 Vec pointVec; 1961 PetscScalar *pV; 1962 PetscSF coarseCellSF = NULL; 1963 const PetscSFNode *coarseCells; 1964 PetscInt numCoarseCells, q, c; 1965 PetscInt *findices, *cindices; 1966 PetscInt numFIndices, numCIndices; 1967 1968 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1969 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1970 /* Get points from the quadrature */ 1971 ierr = VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);CHKERRQ(ierr); 1972 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1973 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1974 for (q = 0; q < Nq; ++q) { 1975 /* Transform point to real space */ 1976 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1977 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1978 } 1979 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1980 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1981 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1982 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 1983 /* Update preallocation info */ 1984 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1985 if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1986 { 1987 PetscHashJKKey key; 1988 PetscHashJKIter missing, iter; 1989 1990 for (i = 0; i < numFIndices; ++i) { 1991 key.j = findices[i]; 1992 if (key.j >= 0) { 1993 /* Get indices for coarse elements */ 1994 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1995 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1996 for (c = 0; c < numCIndices; ++c) { 1997 key.k = cindices[c]; 1998 if (key.k < 0) continue; 1999 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 2000 if (missing) { 2001 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 2002 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 2003 else ++onz[key.j-rStart]; 2004 } 2005 } 2006 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2007 } 2008 } 2009 } 2010 } 2011 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 2012 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 2013 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2014 } 2015 } 2016 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 2017 ierr = MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 2018 ierr = MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2019 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 2020 for (field = 0; field < Nf; ++field) { 2021 PetscObject obj; 2022 PetscClassId id; 2023 PetscQuadrature quad; 2024 PetscReal *Bfine; 2025 const PetscReal *qpoints, *qweights; 2026 PetscInt Nq, Nc, i, d; 2027 2028 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 2029 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2030 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetQuadrature((PetscFE) obj, &quad);CHKERRQ(ierr);ierr = PetscFEGetDefaultTabulation((PetscFE) obj, &Bfine, NULL, NULL);CHKERRQ(ierr);} 2031 else {ierr = PetscFVGetQuadrature((PetscFV) obj, &quad);CHKERRQ(ierr);} 2032 ierr = PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, &qweights);CHKERRQ(ierr); 2033 /* For each fine grid cell */ 2034 for (cell = cStart; cell < cEnd; ++cell) { 2035 Vec pointVec; 2036 PetscScalar *pV; 2037 PetscSF coarseCellSF = NULL; 2038 const PetscSFNode *coarseCells; 2039 PetscInt numCoarseCells, cpdim, q, c, j; 2040 PetscInt *findices, *cindices; 2041 PetscInt numFIndices, numCIndices; 2042 2043 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2044 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 2045 /* Get points from the quadrature */ 2046 ierr = VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);CHKERRQ(ierr); 2047 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 2048 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2049 for (q = 0; q < Nq; ++q) { 2050 /* Transform point to real space */ 2051 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 2052 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 2053 } 2054 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2055 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 2056 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 2057 /* Update matrix */ 2058 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 2059 if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 2060 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2061 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 2062 PetscReal pVReal[3]; 2063 2064 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2065 /* Transform points from real space to coarse reference space */ 2066 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 2067 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 2068 CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x); 2069 2070 if (id == PETSCFE_CLASSID) { 2071 PetscFE fe = (PetscFE) obj; 2072 PetscReal *B; 2073 2074 /* Evaluate coarse basis on contained point */ 2075 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 2076 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 2077 /* Get elemMat entries by multiplying by weight */ 2078 for (i = 0; i < numFIndices; ++i) { 2079 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 2080 for (j = 0; j < cpdim; ++j) { 2081 for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*Bfine[(ccell*numFIndices + i)*Nc + c]*qweights[ccell*Nc + c]*detJ; 2082 } 2083 /* Update interpolator */ 2084 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 2085 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 2086 ierr = MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);CHKERRQ(ierr); 2087 } 2088 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 2089 } else { 2090 cpdim = 1; 2091 for (i = 0; i < numFIndices; ++i) { 2092 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 2093 for (j = 0; j < cpdim; ++j) { 2094 for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*1.0*qweights[ccell*Nc + c]*detJ; 2095 } 2096 /* Update interpolator */ 2097 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 2098 ierr = PetscPrintf(PETSC_COMM_SELF, "Nq: %d %d Nf: %d %d Nc: %d %d\n", ccell, Nq, i, numFIndices, j, numCIndices);CHKERRQ(ierr); 2099 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 2100 ierr = MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);CHKERRQ(ierr); 2101 } 2102 } 2103 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2104 } 2105 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2106 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 2107 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 2108 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2109 } 2110 } 2111 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 2112 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 2113 ierr = PetscFree(elemMat);CHKERRQ(ierr); 2114 ierr = MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2115 ierr = MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2116 PetscFunctionReturn(0); 2117 } 2118 2119 /*@ 2120 DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns 2121 2122 Input Parameters: 2123 + dmc - The coarse mesh 2124 - dmf - The fine mesh 2125 - user - The user context 2126 2127 Output Parameter: 2128 . sc - The mapping 2129 2130 Level: developer 2131 2132 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 2133 @*/ 2134 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 2135 { 2136 PetscDS prob; 2137 PetscFE *feRef; 2138 PetscFV *fvRef; 2139 Vec fv, cv; 2140 IS fis, cis; 2141 PetscSection fsection, fglobalSection, csection, cglobalSection; 2142 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 2143 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m; 2144 PetscErrorCode ierr; 2145 2146 PetscFunctionBegin; 2147 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2148 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 2149 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 2150 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 2151 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 2152 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 2153 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 2154 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 2155 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 2156 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 2157 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 2158 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 2159 for (f = 0; f < Nf; ++f) { 2160 PetscObject obj; 2161 PetscClassId id; 2162 PetscInt fNb = 0, Nc = 0; 2163 2164 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 2165 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2166 if (id == PETSCFE_CLASSID) { 2167 PetscFE fe = (PetscFE) obj; 2168 2169 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 2170 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 2171 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 2172 } else if (id == PETSCFV_CLASSID) { 2173 PetscFV fv = (PetscFV) obj; 2174 PetscDualSpace Q; 2175 2176 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 2177 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 2178 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 2179 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 2180 } 2181 fTotDim += fNb*Nc; 2182 } 2183 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 2184 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 2185 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 2186 PetscFE feC; 2187 PetscFV fvC; 2188 PetscDualSpace QF, QC; 2189 PetscInt NcF, NcC, fpdim, cpdim; 2190 2191 if (feRef[field]) { 2192 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 2193 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 2194 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 2195 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 2196 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 2197 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 2198 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 2199 } else { 2200 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 2201 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 2202 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 2203 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 2204 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 2205 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 2206 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 2207 } 2208 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); 2209 for (c = 0; c < cpdim; ++c) { 2210 PetscQuadrature cfunc; 2211 const PetscReal *cqpoints; 2212 PetscInt NpC; 2213 PetscBool found = PETSC_FALSE; 2214 2215 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 2216 ierr = PetscQuadratureGetData(cfunc, NULL, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 2217 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 2218 for (f = 0; f < fpdim; ++f) { 2219 PetscQuadrature ffunc; 2220 const PetscReal *fqpoints; 2221 PetscReal sum = 0.0; 2222 PetscInt NpF, comp; 2223 2224 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 2225 ierr = PetscQuadratureGetData(ffunc, NULL, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 2226 if (NpC != NpF) continue; 2227 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 2228 if (sum > 1.0e-9) continue; 2229 for (comp = 0; comp < NcC; ++comp) { 2230 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 2231 } 2232 found = PETSC_TRUE; 2233 break; 2234 } 2235 if (!found) { 2236 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 2237 if (fvRef[field]) { 2238 PetscInt comp; 2239 for (comp = 0; comp < NcC; ++comp) { 2240 cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp; 2241 } 2242 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 2243 } 2244 } 2245 offsetC += cpdim*NcC; 2246 offsetF += fpdim*NcF; 2247 } 2248 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 2249 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 2250 2251 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 2252 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 2253 ierr = VecGetOwnershipRange(cv, &startC, &endC);CHKERRQ(ierr); 2254 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 2255 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 2256 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 2257 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 2258 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 2259 for (c = cStart; c < cEnd; ++c) { 2260 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 2261 for (d = 0; d < cTotDim; ++d) { 2262 if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue; 2263 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]]); 2264 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 2265 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 2266 } 2267 } 2268 ierr = PetscFree(cmap);CHKERRQ(ierr); 2269 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 2270 2271 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 2272 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 2273 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 2274 ierr = ISDestroy(&cis);CHKERRQ(ierr); 2275 ierr = ISDestroy(&fis);CHKERRQ(ierr); 2276 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 2277 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 2278 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2279 PetscFunctionReturn(0); 2280 } 2281