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 = NULL; 1153 PetscSection section, sectionAux; 1154 Vec locX, locA; 1155 PetscInt dim, numCells = cEnd - cStart, c, f; 1156 PetscBool 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 PetscFEGeom *cgeomFEM; 1164 DM dmGrad; 1165 PetscQuadrature affineQuad = NULL; 1166 Vec cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL; 1167 PetscFVCellGeom *cgeomFVM; 1168 const PetscScalar *lgrad; 1169 PetscBool isAffine; 1170 DMField coordField; 1171 IS cellIS; 1172 PetscErrorCode ierr; 1173 1174 PetscFunctionBegin; 1175 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1176 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1177 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1178 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1179 /* Determine which discretizations we have */ 1180 for (f = 0; f < Nf; ++f) { 1181 PetscObject obj; 1182 PetscClassId id; 1183 1184 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1185 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1186 if (id == PETSCFV_CLASSID) useFVM = PETSC_TRUE; 1187 } 1188 /* Get local solution with boundary values */ 1189 ierr = DMGetLocalVector(dm, &locX);CHKERRQ(ierr); 1190 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 1191 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);CHKERRQ(ierr); 1192 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);CHKERRQ(ierr); 1193 /* Read DS information */ 1194 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1195 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 1196 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 1197 ierr = ISCreateStride(PETSC_COMM_SELF,numCells,cStart,1,&cellIS);CHKERRQ(ierr); 1198 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 1199 /* Read Auxiliary DS information */ 1200 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1201 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);CHKERRQ(ierr); 1202 if (dmAux) { 1203 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1204 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 1205 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1206 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1207 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 1208 } 1209 /* Allocate data arrays */ 1210 ierr = PetscCalloc1(numCells*totDim, &u);CHKERRQ(ierr); 1211 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1212 /* Read out geometry */ 1213 ierr = DMGetCoordinateField(dm,&coordField);CHKERRQ(ierr); 1214 ierr = DMFieldGetFEInvariance(coordField,cellIS,NULL,&isAffine,NULL);CHKERRQ(ierr); 1215 if (isAffine) { 1216 ierr = DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);CHKERRQ(ierr); 1217 if (affineQuad) { 1218 ierr = DMFieldCreateFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&cgeomFEM);CHKERRQ(ierr); 1219 } 1220 } 1221 if (useFVM) { 1222 PetscFV fv = NULL; 1223 Vec grad; 1224 PetscInt fStart, fEnd; 1225 PetscBool compGrad; 1226 1227 for (f = 0; f < Nf; ++f) { 1228 PetscObject obj; 1229 PetscClassId id; 1230 1231 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1232 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1233 if (id == PETSCFV_CLASSID) {fv = (PetscFV) obj; break;} 1234 } 1235 ierr = PetscFVGetComputeGradients(fv, &compGrad);CHKERRQ(ierr); 1236 ierr = PetscFVSetComputeGradients(fv, PETSC_TRUE);CHKERRQ(ierr); 1237 ierr = DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);CHKERRQ(ierr); 1238 ierr = DMPlexComputeGradientFVM(dm, fv, faceGeometryFVM, cellGeometryFVM, &dmGrad);CHKERRQ(ierr); 1239 ierr = PetscFVSetComputeGradients(fv, compGrad);CHKERRQ(ierr); 1240 ierr = VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1241 /* Reconstruct and limit cell gradients */ 1242 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1243 ierr = DMGetGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1244 ierr = DMPlexReconstructGradients_Internal(dm, fv, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);CHKERRQ(ierr); 1245 /* Communicate gradient values */ 1246 ierr = DMGetLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1247 ierr = DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1248 ierr = DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1249 ierr = DMRestoreGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1250 /* Handle non-essential (e.g. outflow) boundary values */ 1251 ierr = DMPlexInsertBoundaryValues(dm, PETSC_FALSE, locX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);CHKERRQ(ierr); 1252 ierr = VecGetArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1253 } 1254 /* Read out data from inputs */ 1255 for (c = cStart; c < cEnd; ++c) { 1256 PetscScalar *x = NULL; 1257 PetscInt i; 1258 1259 ierr = DMPlexVecGetClosure(dm, section, locX, c, NULL, &x);CHKERRQ(ierr); 1260 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 1261 ierr = DMPlexVecRestoreClosure(dm, section, locX, c, NULL, &x);CHKERRQ(ierr); 1262 if (dmAux) { 1263 ierr = DMPlexVecGetClosure(dmAux, sectionAux, locA, c, NULL, &x);CHKERRQ(ierr); 1264 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 1265 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, locA, c, NULL, &x);CHKERRQ(ierr); 1266 } 1267 } 1268 /* Do integration for each field */ 1269 for (f = 0; f < Nf; ++f) { 1270 PetscObject obj; 1271 PetscClassId id; 1272 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 1273 1274 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1275 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1276 if (id == PETSCFE_CLASSID) { 1277 PetscFE fe = (PetscFE) obj; 1278 PetscQuadrature q; 1279 PetscFEGeom *chunkGeom = NULL; 1280 PetscInt Nq, Nb; 1281 1282 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1283 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1284 ierr = PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 1285 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1286 blockSize = Nb*Nq; 1287 batchSize = numBlocks * blockSize; 1288 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1289 numChunks = numCells / (numBatches*batchSize); 1290 Ne = numChunks*numBatches*batchSize; 1291 Nr = numCells % (numBatches*batchSize); 1292 offset = numCells - Nr; 1293 if (!affineQuad) { 1294 ierr = DMFieldCreateFEGeom(coordField,cellIS,q,PETSC_FALSE,&cgeomFEM);CHKERRQ(ierr); 1295 } 1296 ierr = PetscFEGeomGetChunk(cgeomFEM,0,offset,&chunkGeom);CHKERRQ(ierr); 1297 ierr = PetscFEIntegrate(fe, prob, f, Ne, chunkGeom, u, probAux, a, cintegral);CHKERRQ(ierr); 1298 ierr = PetscFEGeomGetChunk(cgeomFEM,offset,numCells,&chunkGeom);CHKERRQ(ierr); 1299 ierr = PetscFEIntegrate(fe, prob, f, Nr, chunkGeom, &u[offset*totDim], probAux, &a[offset*totDimAux], &cintegral[offset*Nf]);CHKERRQ(ierr); 1300 ierr = PetscFEGeomRestoreChunk(cgeomFEM,offset,numCells,&chunkGeom);CHKERRQ(ierr); 1301 if (!affineQuad) { 1302 ierr = PetscFEGeomDestroy(&cgeomFEM);CHKERRQ(ierr); 1303 } 1304 } else if (id == PETSCFV_CLASSID) { 1305 PetscInt foff; 1306 PetscPointFunc obj_func; 1307 PetscScalar lint; 1308 1309 ierr = PetscDSGetObjective(prob, f, &obj_func);CHKERRQ(ierr); 1310 ierr = PetscDSGetFieldOffset(prob, f, &foff);CHKERRQ(ierr); 1311 if (obj_func) { 1312 for (c = 0; c < numCells; ++c) { 1313 PetscScalar *u_x; 1314 1315 ierr = DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);CHKERRQ(ierr); 1316 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); 1317 cintegral[c*Nf+f] += PetscRealPart(lint)*cgeomFVM[c].volume; 1318 } 1319 } 1320 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 1321 } 1322 /* Cleanup data arrays */ 1323 if (useFVM) { 1324 ierr = VecRestoreArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1325 ierr = VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1326 ierr = DMRestoreLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1327 ierr = VecDestroy(&faceGeometryFVM);CHKERRQ(ierr); 1328 ierr = VecDestroy(&cellGeometryFVM);CHKERRQ(ierr); 1329 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr); 1330 } 1331 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1332 ierr = PetscFree(u);CHKERRQ(ierr); 1333 /* Cleanup */ 1334 if (affineQuad) { 1335 ierr = PetscFEGeomDestroy(&cgeomFEM);CHKERRQ(ierr); 1336 } 1337 ierr = PetscQuadratureDestroy(&affineQuad);CHKERRQ(ierr); 1338 ierr = ISDestroy(&cellIS);CHKERRQ(ierr); 1339 ierr = DMRestoreLocalVector(dm, &locX);CHKERRQ(ierr); 1340 PetscFunctionReturn(0); 1341 } 1342 1343 /*@ 1344 DMPlexComputeIntegralFEM - Form the integral over the domain from the global input X using pointwise functions specified by the user 1345 1346 Input Parameters: 1347 + dm - The mesh 1348 . X - Global input vector 1349 - user - The user context 1350 1351 Output Parameter: 1352 . integral - Integral for each field 1353 1354 Level: developer 1355 1356 .seealso: DMPlexComputeResidualFEM() 1357 @*/ 1358 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscScalar *integral, void *user) 1359 { 1360 DM_Plex *mesh = (DM_Plex *) dm->data; 1361 PetscScalar *cintegral, *lintegral; 1362 PetscInt Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell; 1363 PetscErrorCode ierr; 1364 1365 PetscFunctionBegin; 1366 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1367 PetscValidHeaderSpecific(X, VEC_CLASSID, 2); 1368 PetscValidPointer(integral, 3); 1369 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1370 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 1371 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr); 1372 ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr); 1373 ierr = DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);CHKERRQ(ierr); 1374 cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight]; 1375 /* TODO Introduce a loop over large chunks (right now this is a single chunk) */ 1376 ierr = PetscCalloc2(Nf, &lintegral, (cEnd-cStart)*Nf, &cintegral);CHKERRQ(ierr); 1377 ierr = DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);CHKERRQ(ierr); 1378 /* Sum up values */ 1379 for (cell = cStart; cell < cEnd; ++cell) { 1380 const PetscInt c = cell - cStart; 1381 1382 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);CHKERRQ(ierr);} 1383 for (f = 0; f < Nf; ++f) lintegral[f] += cintegral[c*Nf+f]; 1384 } 1385 ierr = MPIU_Allreduce(lintegral, integral, Nf, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr); 1386 if (mesh->printFEM) { 1387 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Integral:");CHKERRQ(ierr); 1388 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", (double) PetscRealPart(integral[f]));CHKERRQ(ierr);} 1389 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 1390 } 1391 ierr = PetscFree2(lintegral, cintegral);CHKERRQ(ierr); 1392 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1393 PetscFunctionReturn(0); 1394 } 1395 1396 /*@ 1397 DMPlexComputeCellwiseIntegralFEM - Form the vector of cellwise integrals F from the global input X using pointwise functions specified by the user 1398 1399 Input Parameters: 1400 + dm - The mesh 1401 . X - Global input vector 1402 - user - The user context 1403 1404 Output Parameter: 1405 . integral - Cellwise integrals for each field 1406 1407 Level: developer 1408 1409 .seealso: DMPlexComputeResidualFEM() 1410 @*/ 1411 PetscErrorCode DMPlexComputeCellwiseIntegralFEM(DM dm, Vec X, Vec F, void *user) 1412 { 1413 DM_Plex *mesh = (DM_Plex *) dm->data; 1414 DM dmF; 1415 PetscSection sectionF; 1416 PetscScalar *cintegral, *af; 1417 PetscInt Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell; 1418 PetscErrorCode ierr; 1419 1420 PetscFunctionBegin; 1421 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1422 PetscValidHeaderSpecific(X, VEC_CLASSID, 2); 1423 PetscValidHeaderSpecific(F, VEC_CLASSID, 3); 1424 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1425 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 1426 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr); 1427 ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr); 1428 ierr = DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);CHKERRQ(ierr); 1429 cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight]; 1430 /* TODO Introduce a loop over large chunks (right now this is a single chunk) */ 1431 ierr = PetscCalloc1((cEnd-cStart)*Nf, &cintegral);CHKERRQ(ierr); 1432 ierr = DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);CHKERRQ(ierr); 1433 /* Put values in F*/ 1434 ierr = VecGetDM(F, &dmF);CHKERRQ(ierr); 1435 ierr = DMGetDefaultSection(dmF, §ionF);CHKERRQ(ierr); 1436 ierr = VecGetArray(F, &af);CHKERRQ(ierr); 1437 for (cell = cStart; cell < cEnd; ++cell) { 1438 const PetscInt c = cell - cStart; 1439 PetscInt dof, off; 1440 1441 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);CHKERRQ(ierr);} 1442 ierr = PetscSectionGetDof(sectionF, cell, &dof);CHKERRQ(ierr); 1443 ierr = PetscSectionGetOffset(sectionF, cell, &off);CHKERRQ(ierr); 1444 if (dof != Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cell dofs %D != %D", dof, Nf); 1445 for (f = 0; f < Nf; ++f) af[off+f] = cintegral[c*Nf+f]; 1446 } 1447 ierr = VecRestoreArray(F, &af);CHKERRQ(ierr); 1448 ierr = PetscFree(cintegral);CHKERRQ(ierr); 1449 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1450 PetscFunctionReturn(0); 1451 } 1452 1453 /*@ 1454 DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1455 1456 Input Parameters: 1457 + dmf - The fine mesh 1458 . dmc - The coarse mesh 1459 - user - The user context 1460 1461 Output Parameter: 1462 . In - The interpolation matrix 1463 1464 Level: developer 1465 1466 .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1467 @*/ 1468 PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user) 1469 { 1470 DM_Plex *mesh = (DM_Plex *) dmc->data; 1471 const char *name = "Interpolator"; 1472 PetscDS prob; 1473 PetscFE *feRef; 1474 PetscFV *fvRef; 1475 PetscSection fsection, fglobalSection; 1476 PetscSection csection, cglobalSection; 1477 PetscScalar *elemMat; 1478 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c; 1479 PetscInt cTotDim, rTotDim = 0; 1480 PetscErrorCode ierr; 1481 1482 PetscFunctionBegin; 1483 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1484 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1485 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1486 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1487 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1488 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1489 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1490 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1491 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1492 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1493 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr); 1494 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1495 for (f = 0; f < Nf; ++f) { 1496 PetscObject obj; 1497 PetscClassId id; 1498 PetscInt rNb = 0, Nc = 0; 1499 1500 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1501 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1502 if (id == PETSCFE_CLASSID) { 1503 PetscFE fe = (PetscFE) obj; 1504 1505 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1506 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1507 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1508 } else if (id == PETSCFV_CLASSID) { 1509 PetscFV fv = (PetscFV) obj; 1510 PetscDualSpace Q; 1511 1512 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1513 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1514 ierr = PetscDualSpaceGetDimension(Q, &rNb);CHKERRQ(ierr); 1515 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1516 } 1517 rTotDim += rNb; 1518 } 1519 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1520 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 1521 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 1522 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1523 PetscDualSpace Qref; 1524 PetscQuadrature f; 1525 const PetscReal *qpoints, *qweights; 1526 PetscReal *points; 1527 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1528 1529 /* Compose points from all dual basis functionals */ 1530 if (feRef[fieldI]) { 1531 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1532 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 1533 } else { 1534 ierr = PetscFVGetDualSpace(fvRef[fieldI], &Qref);CHKERRQ(ierr); 1535 ierr = PetscFVGetNumComponents(fvRef[fieldI], &Nc);CHKERRQ(ierr); 1536 } 1537 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1538 for (i = 0; i < fpdim; ++i) { 1539 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1540 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1541 npoints += Np; 1542 } 1543 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1544 for (i = 0, k = 0; i < fpdim; ++i) { 1545 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1546 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1547 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1548 } 1549 1550 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1551 PetscObject obj; 1552 PetscClassId id; 1553 PetscReal *B; 1554 PetscInt NcJ = 0, cpdim = 0, j, qNc; 1555 1556 ierr = PetscDSGetDiscretization(prob, fieldJ, &obj);CHKERRQ(ierr); 1557 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1558 if (id == PETSCFE_CLASSID) { 1559 PetscFE fe = (PetscFE) obj; 1560 1561 /* Evaluate basis at points */ 1562 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 1563 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1564 /* For now, fields only interpolate themselves */ 1565 if (fieldI == fieldJ) { 1566 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); 1567 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1568 for (i = 0, k = 0; i < fpdim; ++i) { 1569 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1570 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);CHKERRQ(ierr); 1571 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); 1572 for (p = 0; p < Np; ++p, ++k) { 1573 for (j = 0; j < cpdim; ++j) { 1574 /* 1575 cTotDim: Total columns in element interpolation matrix, sum of number of dual basis functionals in each field 1576 offsetI, offsetJ: Offsets into the larger element interpolation matrix for different fields 1577 fpdim, i, cpdim, j: Dofs for fine and coarse grids, correspond to dual space basis functionals 1578 qNC, Nc, Ncj, c: Number of components in this field 1579 Np, p: Number of quad points in the fine grid functional i 1580 k: i*Np + p, overall point number for the interpolation 1581 */ 1582 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p*qNc+c]; 1583 } 1584 } 1585 } 1586 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1587 } 1588 } else if (id == PETSCFV_CLASSID) { 1589 PetscFV fv = (PetscFV) obj; 1590 1591 /* Evaluate constant function at points */ 1592 ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr); 1593 cpdim = 1; 1594 /* For now, fields only interpolate themselves */ 1595 if (fieldI == fieldJ) { 1596 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); 1597 for (i = 0, k = 0; i < fpdim; ++i) { 1598 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1599 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);CHKERRQ(ierr); 1600 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); 1601 for (p = 0; p < Np; ++p, ++k) { 1602 for (j = 0; j < cpdim; ++j) { 1603 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p*qNc+c]; 1604 } 1605 } 1606 } 1607 } 1608 } 1609 offsetJ += cpdim; 1610 } 1611 offsetI += fpdim; 1612 ierr = PetscFree(points);CHKERRQ(ierr); 1613 } 1614 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1615 /* Preallocate matrix */ 1616 { 1617 Mat preallocator; 1618 PetscScalar *vals; 1619 PetscInt *cellCIndices, *cellFIndices; 1620 PetscInt locRows, locCols, cell; 1621 1622 ierr = MatGetLocalSize(In, &locRows, &locCols);CHKERRQ(ierr); 1623 ierr = MatCreate(PetscObjectComm((PetscObject) In), &preallocator);CHKERRQ(ierr); 1624 ierr = MatSetType(preallocator, MATPREALLOCATOR);CHKERRQ(ierr); 1625 ierr = MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1626 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 1627 ierr = PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1628 for (cell = cStart; cell < cEnd; ++cell) { 1629 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1630 ierr = MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);CHKERRQ(ierr); 1631 } 1632 ierr = PetscFree3(vals,cellCIndices,cellFIndices);CHKERRQ(ierr); 1633 ierr = MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1634 ierr = MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1635 ierr = MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);CHKERRQ(ierr); 1636 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 1637 } 1638 /* Fill matrix */ 1639 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1640 for (c = cStart; c < cEnd; ++c) { 1641 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1642 } 1643 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1644 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1645 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1646 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1647 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1648 if (mesh->printFEM) { 1649 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1650 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1651 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1652 } 1653 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1654 PetscFunctionReturn(0); 1655 } 1656 1657 PetscErrorCode DMPlexComputeMassMatrixNested(DM dmc, DM dmf, Mat mass, void *user) 1658 { 1659 SETERRQ(PetscObjectComm((PetscObject) dmc), PETSC_ERR_SUP, "Laziness"); 1660 } 1661 1662 /*@ 1663 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 1664 1665 Input Parameters: 1666 + dmf - The fine mesh 1667 . dmc - The coarse mesh 1668 - user - The user context 1669 1670 Output Parameter: 1671 . In - The interpolation matrix 1672 1673 Level: developer 1674 1675 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1676 @*/ 1677 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 1678 { 1679 DM_Plex *mesh = (DM_Plex *) dmf->data; 1680 const char *name = "Interpolator"; 1681 PetscDS prob; 1682 PetscSection fsection, csection, globalFSection, globalCSection; 1683 PetscHashJK ht; 1684 PetscLayout rLayout; 1685 PetscInt *dnz, *onz; 1686 PetscInt locRows, rStart, rEnd; 1687 PetscReal *x, *v0, *J, *invJ, detJ; 1688 PetscReal *v0c, *Jc, *invJc, detJc; 1689 PetscScalar *elemMat; 1690 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1691 PetscErrorCode ierr; 1692 1693 PetscFunctionBegin; 1694 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1695 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1696 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1697 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1698 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1699 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1700 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1701 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1702 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1703 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1704 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1705 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1706 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1707 ierr = PetscMalloc1(totDim, &elemMat);CHKERRQ(ierr); 1708 1709 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1710 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1711 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1712 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1713 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1714 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1715 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1716 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1717 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1718 for (field = 0; field < Nf; ++field) { 1719 PetscObject obj; 1720 PetscClassId id; 1721 PetscDualSpace Q = NULL; 1722 PetscQuadrature f; 1723 const PetscReal *qpoints; 1724 PetscInt Nc, Np, fpdim, i, d; 1725 1726 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1727 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1728 if (id == PETSCFE_CLASSID) { 1729 PetscFE fe = (PetscFE) obj; 1730 1731 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1732 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1733 } else if (id == PETSCFV_CLASSID) { 1734 PetscFV fv = (PetscFV) obj; 1735 1736 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1737 Nc = 1; 1738 } 1739 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1740 /* For each fine grid cell */ 1741 for (cell = cStart; cell < cEnd; ++cell) { 1742 PetscInt *findices, *cindices; 1743 PetscInt numFIndices, numCIndices; 1744 1745 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1746 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1747 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim); 1748 for (i = 0; i < fpdim; ++i) { 1749 Vec pointVec; 1750 PetscScalar *pV; 1751 PetscSF coarseCellSF = NULL; 1752 const PetscSFNode *coarseCells; 1753 PetscInt numCoarseCells, q, c; 1754 1755 /* Get points from the dual basis functional quadrature */ 1756 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1757 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1758 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1759 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1760 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1761 for (q = 0; q < Np; ++q) { 1762 const PetscReal xi0[3] = {-1., -1., -1.}; 1763 1764 /* Transform point to real space */ 1765 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x); 1766 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1767 } 1768 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1769 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1770 /* OPT: Pack all quad points from fine cell */ 1771 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1772 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 1773 /* Update preallocation info */ 1774 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1775 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1776 { 1777 PetscHashJKKey key; 1778 PetscHashJKIter missing, iter; 1779 1780 key.j = findices[i]; 1781 if (key.j >= 0) { 1782 /* Get indices for coarse elements */ 1783 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1784 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1785 for (c = 0; c < numCIndices; ++c) { 1786 key.k = cindices[c]; 1787 if (key.k < 0) continue; 1788 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1789 if (missing) { 1790 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1791 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1792 else ++onz[key.j-rStart]; 1793 } 1794 } 1795 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1796 } 1797 } 1798 } 1799 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1800 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1801 } 1802 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1803 } 1804 } 1805 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1806 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1807 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1808 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1809 for (field = 0; field < Nf; ++field) { 1810 PetscObject obj; 1811 PetscClassId id; 1812 PetscDualSpace Q = NULL; 1813 PetscQuadrature f; 1814 const PetscReal *qpoints, *qweights; 1815 PetscInt Nc, qNc, Np, fpdim, i, d; 1816 1817 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1818 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1819 if (id == PETSCFE_CLASSID) { 1820 PetscFE fe = (PetscFE) obj; 1821 1822 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1823 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1824 } else if (id == PETSCFV_CLASSID) { 1825 PetscFV fv = (PetscFV) obj; 1826 1827 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1828 Nc = 1; 1829 } else SETERRQ1(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONG,"Unknown discretization type for field %d",field); 1830 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1831 /* For each fine grid cell */ 1832 for (cell = cStart; cell < cEnd; ++cell) { 1833 PetscInt *findices, *cindices; 1834 PetscInt numFIndices, numCIndices; 1835 1836 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1837 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1838 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim); 1839 for (i = 0; i < fpdim; ++i) { 1840 Vec pointVec; 1841 PetscScalar *pV; 1842 PetscSF coarseCellSF = NULL; 1843 const PetscSFNode *coarseCells; 1844 PetscInt numCoarseCells, cpdim, q, c, j; 1845 1846 /* Get points from the dual basis functional quadrature */ 1847 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1848 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);CHKERRQ(ierr); 1849 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); 1850 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1851 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1852 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1853 for (q = 0; q < Np; ++q) { 1854 const PetscReal xi0[3] = {-1., -1., -1.}; 1855 1856 /* Transform point to real space */ 1857 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x); 1858 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1859 } 1860 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1861 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1862 /* OPT: Read this out from preallocation information */ 1863 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1864 /* Update preallocation info */ 1865 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1866 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1867 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1868 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1869 PetscReal pVReal[3]; 1870 const PetscReal xi0[3] = {-1., -1., -1.}; 1871 1872 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1873 /* Transform points from real space to coarse reference space */ 1874 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 1875 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 1876 CoordinatesRealToRef(dim, dim, xi0, v0c, invJc, pVReal, x); 1877 1878 if (id == PETSCFE_CLASSID) { 1879 PetscFE fe = (PetscFE) obj; 1880 PetscReal *B; 1881 1882 /* Evaluate coarse basis on contained point */ 1883 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1884 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 1885 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 1886 /* Get elemMat entries by multiplying by weight */ 1887 for (j = 0; j < cpdim; ++j) { 1888 for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*qweights[ccell*qNc + c]; 1889 } 1890 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1891 } else { 1892 cpdim = 1; 1893 for (j = 0; j < cpdim; ++j) { 1894 for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c]; 1895 } 1896 } 1897 /* Update interpolator */ 1898 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 1899 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 1900 ierr = MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1901 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1902 } 1903 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1904 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1905 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1906 } 1907 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1908 } 1909 } 1910 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 1911 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 1912 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1913 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1914 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1915 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1916 PetscFunctionReturn(0); 1917 } 1918 1919 /*@ 1920 DMPlexComputeMassMatrixGeneral - Form the local portion of the mass matrix M from the coarse DM to a non-nested fine DM. 1921 1922 Input Parameters: 1923 + dmf - The fine mesh 1924 . dmc - The coarse mesh 1925 - user - The user context 1926 1927 Output Parameter: 1928 . mass - The mass matrix 1929 1930 Level: developer 1931 1932 .seealso: DMPlexComputeMassMatrixNested(), DMPlexComputeInterpolatorNested(), DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1933 @*/ 1934 PetscErrorCode DMPlexComputeMassMatrixGeneral(DM dmc, DM dmf, Mat mass, void *user) 1935 { 1936 DM_Plex *mesh = (DM_Plex *) dmf->data; 1937 const char *name = "Mass Matrix"; 1938 PetscDS prob; 1939 PetscSection fsection, csection, globalFSection, globalCSection; 1940 PetscHashJK ht; 1941 PetscLayout rLayout; 1942 PetscInt *dnz, *onz; 1943 PetscInt locRows, rStart, rEnd; 1944 PetscReal *x, *v0, *J, *invJ, detJ; 1945 PetscReal *v0c, *Jc, *invJc, detJc; 1946 PetscScalar *elemMat; 1947 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1948 PetscErrorCode ierr; 1949 1950 PetscFunctionBegin; 1951 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1952 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1953 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1954 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1955 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1956 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1957 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1958 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1959 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1960 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1961 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1962 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1963 ierr = PetscMalloc1(totDim, &elemMat);CHKERRQ(ierr); 1964 1965 ierr = MatGetLocalSize(mass, &locRows, NULL);CHKERRQ(ierr); 1966 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) mass), &rLayout);CHKERRQ(ierr); 1967 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1968 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1969 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1970 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1971 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1972 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1973 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1974 for (field = 0; field < Nf; ++field) { 1975 PetscObject obj; 1976 PetscClassId id; 1977 PetscQuadrature quad; 1978 const PetscReal *qpoints; 1979 PetscInt Nq, Nc, i, d; 1980 1981 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1982 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1983 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetQuadrature((PetscFE) obj, &quad);CHKERRQ(ierr);} 1984 else {ierr = PetscFVGetQuadrature((PetscFV) obj, &quad);CHKERRQ(ierr);} 1985 ierr = PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, NULL);CHKERRQ(ierr); 1986 /* For each fine grid cell */ 1987 for (cell = cStart; cell < cEnd; ++cell) { 1988 Vec pointVec; 1989 PetscScalar *pV; 1990 PetscSF coarseCellSF = NULL; 1991 const PetscSFNode *coarseCells; 1992 PetscInt numCoarseCells, q, c; 1993 PetscInt *findices, *cindices; 1994 PetscInt numFIndices, numCIndices; 1995 1996 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1997 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1998 /* Get points from the quadrature */ 1999 ierr = VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);CHKERRQ(ierr); 2000 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 2001 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2002 for (q = 0; q < Nq; ++q) { 2003 const PetscReal xi0[3] = {-1., -1., -1.}; 2004 2005 /* Transform point to real space */ 2006 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x); 2007 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 2008 } 2009 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2010 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 2011 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 2012 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 2013 /* Update preallocation info */ 2014 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 2015 if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 2016 { 2017 PetscHashJKKey key; 2018 PetscHashJKIter missing, iter; 2019 2020 for (i = 0; i < numFIndices; ++i) { 2021 key.j = findices[i]; 2022 if (key.j >= 0) { 2023 /* Get indices for coarse elements */ 2024 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 2025 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2026 for (c = 0; c < numCIndices; ++c) { 2027 key.k = cindices[c]; 2028 if (key.k < 0) continue; 2029 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 2030 if (missing) { 2031 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 2032 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 2033 else ++onz[key.j-rStart]; 2034 } 2035 } 2036 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2037 } 2038 } 2039 } 2040 } 2041 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 2042 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 2043 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2044 } 2045 } 2046 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 2047 ierr = MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 2048 ierr = MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2049 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 2050 for (field = 0; field < Nf; ++field) { 2051 PetscObject obj; 2052 PetscClassId id; 2053 PetscQuadrature quad; 2054 PetscReal *Bfine; 2055 const PetscReal *qpoints, *qweights; 2056 PetscInt Nq, Nc, i, d; 2057 2058 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 2059 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2060 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetQuadrature((PetscFE) obj, &quad);CHKERRQ(ierr);ierr = PetscFEGetDefaultTabulation((PetscFE) obj, &Bfine, NULL, NULL);CHKERRQ(ierr);} 2061 else {ierr = PetscFVGetQuadrature((PetscFV) obj, &quad);CHKERRQ(ierr);} 2062 ierr = PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, &qweights);CHKERRQ(ierr); 2063 /* For each fine grid cell */ 2064 for (cell = cStart; cell < cEnd; ++cell) { 2065 Vec pointVec; 2066 PetscScalar *pV; 2067 PetscSF coarseCellSF = NULL; 2068 const PetscSFNode *coarseCells; 2069 PetscInt numCoarseCells, cpdim, q, c, j; 2070 PetscInt *findices, *cindices; 2071 PetscInt numFIndices, numCIndices; 2072 2073 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2074 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 2075 /* Get points from the quadrature */ 2076 ierr = VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);CHKERRQ(ierr); 2077 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 2078 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2079 for (q = 0; q < Nq; ++q) { 2080 const PetscReal xi0[3] = {-1., -1., -1.}; 2081 2082 /* Transform point to real space */ 2083 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x); 2084 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 2085 } 2086 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2087 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 2088 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 2089 /* Update matrix */ 2090 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 2091 if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 2092 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2093 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 2094 PetscReal pVReal[3]; 2095 const PetscReal xi0[3] = {-1., -1., -1.}; 2096 2097 2098 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2099 /* Transform points from real space to coarse reference space */ 2100 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 2101 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 2102 CoordinatesRealToRef(dim, dim, xi0, v0c, invJc, pVReal, x); 2103 2104 if (id == PETSCFE_CLASSID) { 2105 PetscFE fe = (PetscFE) obj; 2106 PetscReal *B; 2107 2108 /* Evaluate coarse basis on contained point */ 2109 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 2110 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 2111 /* Get elemMat entries by multiplying by weight */ 2112 for (i = 0; i < numFIndices; ++i) { 2113 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 2114 for (j = 0; j < cpdim; ++j) { 2115 for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*Bfine[(ccell*numFIndices + i)*Nc + c]*qweights[ccell*Nc + c]*detJ; 2116 } 2117 /* Update interpolator */ 2118 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 2119 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 2120 ierr = MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);CHKERRQ(ierr); 2121 } 2122 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 2123 } else { 2124 cpdim = 1; 2125 for (i = 0; i < numFIndices; ++i) { 2126 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 2127 for (j = 0; j < cpdim; ++j) { 2128 for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*1.0*qweights[ccell*Nc + c]*detJ; 2129 } 2130 /* Update interpolator */ 2131 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 2132 ierr = PetscPrintf(PETSC_COMM_SELF, "Nq: %d %d Nf: %d %d Nc: %d %d\n", ccell, Nq, i, numFIndices, j, numCIndices);CHKERRQ(ierr); 2133 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 2134 ierr = MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);CHKERRQ(ierr); 2135 } 2136 } 2137 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2138 } 2139 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2140 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 2141 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 2142 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2143 } 2144 } 2145 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 2146 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 2147 ierr = PetscFree(elemMat);CHKERRQ(ierr); 2148 ierr = MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2149 ierr = MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2150 PetscFunctionReturn(0); 2151 } 2152 2153 /*@ 2154 DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns 2155 2156 Input Parameters: 2157 + dmc - The coarse mesh 2158 - dmf - The fine mesh 2159 - user - The user context 2160 2161 Output Parameter: 2162 . sc - The mapping 2163 2164 Level: developer 2165 2166 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 2167 @*/ 2168 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 2169 { 2170 PetscDS prob; 2171 PetscFE *feRef; 2172 PetscFV *fvRef; 2173 Vec fv, cv; 2174 IS fis, cis; 2175 PetscSection fsection, fglobalSection, csection, cglobalSection; 2176 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 2177 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m; 2178 PetscErrorCode ierr; 2179 2180 PetscFunctionBegin; 2181 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2182 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 2183 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 2184 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 2185 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 2186 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 2187 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 2188 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 2189 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 2190 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 2191 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 2192 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 2193 for (f = 0; f < Nf; ++f) { 2194 PetscObject obj; 2195 PetscClassId id; 2196 PetscInt fNb = 0, Nc = 0; 2197 2198 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 2199 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2200 if (id == PETSCFE_CLASSID) { 2201 PetscFE fe = (PetscFE) obj; 2202 2203 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 2204 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 2205 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 2206 } else if (id == PETSCFV_CLASSID) { 2207 PetscFV fv = (PetscFV) obj; 2208 PetscDualSpace Q; 2209 2210 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 2211 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 2212 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 2213 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 2214 } 2215 fTotDim += fNb; 2216 } 2217 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 2218 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 2219 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 2220 PetscFE feC; 2221 PetscFV fvC; 2222 PetscDualSpace QF, QC; 2223 PetscInt order = -1, NcF, NcC, fpdim, cpdim; 2224 2225 if (feRef[field]) { 2226 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 2227 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 2228 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 2229 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 2230 ierr = PetscDualSpaceGetOrder(QF, &order);CHKERRQ(ierr); 2231 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 2232 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 2233 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 2234 } else { 2235 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 2236 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 2237 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 2238 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 2239 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 2240 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 2241 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 2242 } 2243 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); 2244 for (c = 0; c < cpdim; ++c) { 2245 PetscQuadrature cfunc; 2246 const PetscReal *cqpoints, *cqweights; 2247 PetscInt NqcC, NpC; 2248 PetscBool found = PETSC_FALSE; 2249 2250 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 2251 ierr = PetscQuadratureGetData(cfunc, NULL, &NqcC, &NpC, &cqpoints, &cqweights);CHKERRQ(ierr); 2252 if (NqcC != NcC) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of quadrature components %D must match number of field components", NqcC, NcC); 2253 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 2254 for (f = 0; f < fpdim; ++f) { 2255 PetscQuadrature ffunc; 2256 const PetscReal *fqpoints, *fqweights; 2257 PetscReal sum = 0.0; 2258 PetscInt NqcF, NpF; 2259 2260 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 2261 ierr = PetscQuadratureGetData(ffunc, NULL, &NqcF, &NpF, &fqpoints, &fqweights);CHKERRQ(ierr); 2262 if (NqcF != NcF) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of quadrature components %D must match number of field components", NqcF, NcF); 2263 if (NpC != NpF) continue; 2264 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 2265 if (sum > 1.0e-9) continue; 2266 for (d = 0; d < NcC; ++d) sum += PetscAbsReal(cqweights[d]*fqweights[d]); 2267 if (sum < 1.0e-9) continue; 2268 cmap[offsetC+c] = offsetF+f; 2269 found = PETSC_TRUE; 2270 break; 2271 } 2272 if (!found) { 2273 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 2274 if (fvRef[field] || (feRef[field] && order == 0)) { 2275 cmap[offsetC+c] = offsetF+0; 2276 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 2277 } 2278 } 2279 offsetC += cpdim; 2280 offsetF += fpdim; 2281 } 2282 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 2283 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 2284 2285 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 2286 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 2287 ierr = VecGetOwnershipRange(cv, &startC, &endC);CHKERRQ(ierr); 2288 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 2289 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 2290 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 2291 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 2292 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 2293 for (c = cStart; c < cEnd; ++c) { 2294 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 2295 for (d = 0; d < cTotDim; ++d) { 2296 if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue; 2297 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]]); 2298 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 2299 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 2300 } 2301 } 2302 ierr = PetscFree(cmap);CHKERRQ(ierr); 2303 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 2304 2305 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 2306 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 2307 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 2308 ierr = ISDestroy(&cis);CHKERRQ(ierr); 2309 ierr = ISDestroy(&fis);CHKERRQ(ierr); 2310 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 2311 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 2312 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2313 PetscFunctionReturn(0); 2314 } 2315