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 /* NOTE: This is not quite right, unless fpdim == number of fine grid functional quad points */ 1575 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p*qNc+c]; 1576 } 1577 } 1578 } 1579 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1580 } 1581 } else if (id == PETSCFV_CLASSID) { 1582 PetscFV fv = (PetscFV) obj; 1583 1584 /* Evaluate constant function at points */ 1585 ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr); 1586 cpdim = 1; 1587 /* For now, fields only interpolate themselves */ 1588 if (fieldI == fieldJ) { 1589 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); 1590 for (i = 0, k = 0; i < fpdim; ++i) { 1591 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1592 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);CHKERRQ(ierr); 1593 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); 1594 for (p = 0; p < Np; ++p, ++k) { 1595 for (j = 0; j < cpdim; ++j) { 1596 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p*qNc+c]; 1597 } 1598 } 1599 } 1600 } 1601 } 1602 offsetJ += cpdim*NcJ; 1603 } 1604 offsetI += fpdim*Nc; 1605 ierr = PetscFree(points);CHKERRQ(ierr); 1606 } 1607 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1608 /* Preallocate matrix */ 1609 { 1610 Mat preallocator; 1611 PetscScalar *vals; 1612 PetscInt *cellCIndices, *cellFIndices; 1613 PetscInt locRows, locCols, cell; 1614 1615 ierr = MatGetLocalSize(In, &locRows, &locCols);CHKERRQ(ierr); 1616 ierr = MatCreate(PetscObjectComm((PetscObject) In), &preallocator);CHKERRQ(ierr); 1617 ierr = MatSetType(preallocator, MATPREALLOCATOR);CHKERRQ(ierr); 1618 ierr = MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1619 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 1620 ierr = PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1621 for (cell = cStart; cell < cEnd; ++cell) { 1622 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1623 ierr = MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);CHKERRQ(ierr); 1624 } 1625 ierr = PetscFree3(vals,cellCIndices,cellFIndices);CHKERRQ(ierr); 1626 ierr = MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1627 ierr = MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1628 ierr = MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);CHKERRQ(ierr); 1629 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 1630 } 1631 /* Fill matrix */ 1632 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1633 for (c = cStart; c < cEnd; ++c) { 1634 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1635 } 1636 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1637 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1638 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1639 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1640 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1641 if (mesh->printFEM) { 1642 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1643 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1644 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1645 } 1646 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1647 PetscFunctionReturn(0); 1648 } 1649 1650 PetscErrorCode DMPlexComputeMassMatrixNested(DM dmc, DM dmf, Mat mass, void *user) 1651 { 1652 SETERRQ(PetscObjectComm((PetscObject) dmc), PETSC_ERR_SUP, "Laziness"); 1653 } 1654 1655 /*@ 1656 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 1657 1658 Input Parameters: 1659 + dmf - The fine mesh 1660 . dmc - The coarse mesh 1661 - user - The user context 1662 1663 Output Parameter: 1664 . In - The interpolation matrix 1665 1666 Level: developer 1667 1668 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1669 @*/ 1670 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 1671 { 1672 DM_Plex *mesh = (DM_Plex *) dmf->data; 1673 const char *name = "Interpolator"; 1674 PetscDS prob; 1675 PetscSection fsection, csection, globalFSection, globalCSection; 1676 PetscHashJK ht; 1677 PetscLayout rLayout; 1678 PetscInt *dnz, *onz; 1679 PetscInt locRows, rStart, rEnd; 1680 PetscReal *x, *v0, *J, *invJ, detJ; 1681 PetscReal *v0c, *Jc, *invJc, detJc; 1682 PetscScalar *elemMat; 1683 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1684 PetscErrorCode ierr; 1685 1686 PetscFunctionBegin; 1687 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1688 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1689 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1690 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1691 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1692 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1693 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1694 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1695 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1696 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1697 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1698 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1699 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1700 ierr = PetscMalloc1(totDim, &elemMat);CHKERRQ(ierr); 1701 1702 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1703 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1704 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1705 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1706 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1707 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1708 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1709 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1710 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1711 for (field = 0; field < Nf; ++field) { 1712 PetscObject obj; 1713 PetscClassId id; 1714 PetscDualSpace Q = NULL; 1715 PetscQuadrature f; 1716 const PetscReal *qpoints; 1717 PetscInt Nc, Np, fpdim, i, d; 1718 1719 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1720 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1721 if (id == PETSCFE_CLASSID) { 1722 PetscFE fe = (PetscFE) obj; 1723 1724 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1725 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1726 } else if (id == PETSCFV_CLASSID) { 1727 PetscFV fv = (PetscFV) obj; 1728 1729 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1730 Nc = 1; 1731 } 1732 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1733 /* For each fine grid cell */ 1734 for (cell = cStart; cell < cEnd; ++cell) { 1735 PetscInt *findices, *cindices; 1736 PetscInt numFIndices, numCIndices; 1737 1738 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1739 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1740 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim); 1741 for (i = 0; i < fpdim; ++i) { 1742 Vec pointVec; 1743 PetscScalar *pV; 1744 PetscSF coarseCellSF = NULL; 1745 const PetscSFNode *coarseCells; 1746 PetscInt numCoarseCells, q, c; 1747 1748 /* Get points from the dual basis functional quadrature */ 1749 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1750 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1751 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1752 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1753 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1754 for (q = 0; q < Np; ++q) { 1755 const PetscReal xi0[3] = {-1., -1., -1.}; 1756 1757 /* Transform point to real space */ 1758 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x); 1759 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1760 } 1761 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1762 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1763 /* OPT: Pack all quad points from fine cell */ 1764 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1765 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 1766 /* Update preallocation info */ 1767 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1768 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1769 { 1770 PetscHashJKKey key; 1771 PetscHashJKIter missing, iter; 1772 1773 key.j = findices[i]; 1774 if (key.j >= 0) { 1775 /* Get indices for coarse elements */ 1776 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1777 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1778 for (c = 0; c < numCIndices; ++c) { 1779 key.k = cindices[c]; 1780 if (key.k < 0) continue; 1781 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1782 if (missing) { 1783 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1784 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1785 else ++onz[key.j-rStart]; 1786 } 1787 } 1788 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1789 } 1790 } 1791 } 1792 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1793 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1794 } 1795 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1796 } 1797 } 1798 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1799 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1800 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1801 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1802 for (field = 0; field < Nf; ++field) { 1803 PetscObject obj; 1804 PetscClassId id; 1805 PetscDualSpace Q = NULL; 1806 PetscQuadrature f; 1807 const PetscReal *qpoints, *qweights; 1808 PetscInt Nc, qNc, Np, fpdim, i, d; 1809 1810 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1811 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1812 if (id == PETSCFE_CLASSID) { 1813 PetscFE fe = (PetscFE) obj; 1814 1815 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1816 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1817 } else if (id == PETSCFV_CLASSID) { 1818 PetscFV fv = (PetscFV) obj; 1819 1820 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1821 Nc = 1; 1822 } 1823 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1824 /* For each fine grid cell */ 1825 for (cell = cStart; cell < cEnd; ++cell) { 1826 PetscInt *findices, *cindices; 1827 PetscInt numFIndices, numCIndices; 1828 1829 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1830 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1831 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim); 1832 for (i = 0; i < fpdim; ++i) { 1833 Vec pointVec; 1834 PetscScalar *pV; 1835 PetscSF coarseCellSF = NULL; 1836 const PetscSFNode *coarseCells; 1837 PetscInt numCoarseCells, cpdim, q, c, j; 1838 1839 /* Get points from the dual basis functional quadrature */ 1840 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1841 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);CHKERRQ(ierr); 1842 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); 1843 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1844 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1845 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1846 for (q = 0; q < Np; ++q) { 1847 const PetscReal xi0[3] = {-1., -1., -1.}; 1848 1849 /* Transform point to real space */ 1850 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x); 1851 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1852 } 1853 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1854 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1855 /* OPT: Read this out from preallocation information */ 1856 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1857 /* Update preallocation info */ 1858 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1859 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1860 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1861 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1862 PetscReal pVReal[3]; 1863 const PetscReal xi0[3] = {-1., -1., -1.}; 1864 1865 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1866 /* Transform points from real space to coarse reference space */ 1867 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 1868 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 1869 CoordinatesRealToRef(dim, dim, xi0, v0c, invJc, pVReal, x); 1870 1871 if (id == PETSCFE_CLASSID) { 1872 PetscFE fe = (PetscFE) obj; 1873 PetscReal *B; 1874 1875 /* Evaluate coarse basis on contained point */ 1876 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1877 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 1878 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 1879 /* Get elemMat entries by multiplying by weight */ 1880 for (j = 0; j < cpdim; ++j) { 1881 for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*qweights[ccell*qNc + c]; 1882 } 1883 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1884 } else { 1885 cpdim = 1; 1886 for (j = 0; j < cpdim; ++j) { 1887 for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c]; 1888 } 1889 } 1890 /* Update interpolator */ 1891 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 1892 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 1893 ierr = MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1894 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1895 } 1896 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1897 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1898 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1899 } 1900 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1901 } 1902 } 1903 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 1904 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 1905 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1906 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1907 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1908 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1909 PetscFunctionReturn(0); 1910 } 1911 1912 /*@ 1913 DMPlexComputeMassMatrixGeneral - Form the local portion of the mass matrix M from the coarse DM to a non-nested fine DM. 1914 1915 Input Parameters: 1916 + dmf - The fine mesh 1917 . dmc - The coarse mesh 1918 - user - The user context 1919 1920 Output Parameter: 1921 . mass - The mass matrix 1922 1923 Level: developer 1924 1925 .seealso: DMPlexComputeMassMatrixNested(), DMPlexComputeInterpolatorNested(), DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1926 @*/ 1927 PetscErrorCode DMPlexComputeMassMatrixGeneral(DM dmc, DM dmf, Mat mass, void *user) 1928 { 1929 DM_Plex *mesh = (DM_Plex *) dmf->data; 1930 const char *name = "Mass Matrix"; 1931 PetscDS prob; 1932 PetscSection fsection, csection, globalFSection, globalCSection; 1933 PetscHashJK ht; 1934 PetscLayout rLayout; 1935 PetscInt *dnz, *onz; 1936 PetscInt locRows, rStart, rEnd; 1937 PetscReal *x, *v0, *J, *invJ, detJ; 1938 PetscReal *v0c, *Jc, *invJc, detJc; 1939 PetscScalar *elemMat; 1940 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1941 PetscErrorCode ierr; 1942 1943 PetscFunctionBegin; 1944 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1945 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1946 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1947 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1948 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1949 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1950 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1951 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1952 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1953 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1954 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1955 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1956 ierr = PetscMalloc1(totDim, &elemMat);CHKERRQ(ierr); 1957 1958 ierr = MatGetLocalSize(mass, &locRows, NULL);CHKERRQ(ierr); 1959 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) mass), &rLayout);CHKERRQ(ierr); 1960 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1961 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1962 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1963 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1964 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1965 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1966 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1967 for (field = 0; field < Nf; ++field) { 1968 PetscObject obj; 1969 PetscClassId id; 1970 PetscQuadrature quad; 1971 const PetscReal *qpoints; 1972 PetscInt Nq, Nc, i, d; 1973 1974 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1975 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1976 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetQuadrature((PetscFE) obj, &quad);CHKERRQ(ierr);} 1977 else {ierr = PetscFVGetQuadrature((PetscFV) obj, &quad);CHKERRQ(ierr);} 1978 ierr = PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, NULL);CHKERRQ(ierr); 1979 /* For each fine grid cell */ 1980 for (cell = cStart; cell < cEnd; ++cell) { 1981 Vec pointVec; 1982 PetscScalar *pV; 1983 PetscSF coarseCellSF = NULL; 1984 const PetscSFNode *coarseCells; 1985 PetscInt numCoarseCells, q, c; 1986 PetscInt *findices, *cindices; 1987 PetscInt numFIndices, numCIndices; 1988 1989 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1990 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1991 /* Get points from the quadrature */ 1992 ierr = VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);CHKERRQ(ierr); 1993 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1994 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1995 for (q = 0; q < Nq; ++q) { 1996 const PetscReal xi0[3] = {-1., -1., -1.}; 1997 1998 /* Transform point to real space */ 1999 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x); 2000 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 2001 } 2002 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2003 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 2004 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 2005 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 2006 /* Update preallocation info */ 2007 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 2008 if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 2009 { 2010 PetscHashJKKey key; 2011 PetscHashJKIter missing, iter; 2012 2013 for (i = 0; i < numFIndices; ++i) { 2014 key.j = findices[i]; 2015 if (key.j >= 0) { 2016 /* Get indices for coarse elements */ 2017 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 2018 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2019 for (c = 0; c < numCIndices; ++c) { 2020 key.k = cindices[c]; 2021 if (key.k < 0) continue; 2022 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 2023 if (missing) { 2024 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 2025 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 2026 else ++onz[key.j-rStart]; 2027 } 2028 } 2029 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2030 } 2031 } 2032 } 2033 } 2034 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 2035 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 2036 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2037 } 2038 } 2039 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 2040 ierr = MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 2041 ierr = MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2042 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 2043 for (field = 0; field < Nf; ++field) { 2044 PetscObject obj; 2045 PetscClassId id; 2046 PetscQuadrature quad; 2047 PetscReal *Bfine; 2048 const PetscReal *qpoints, *qweights; 2049 PetscInt Nq, Nc, i, d; 2050 2051 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 2052 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2053 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetQuadrature((PetscFE) obj, &quad);CHKERRQ(ierr);ierr = PetscFEGetDefaultTabulation((PetscFE) obj, &Bfine, NULL, NULL);CHKERRQ(ierr);} 2054 else {ierr = PetscFVGetQuadrature((PetscFV) obj, &quad);CHKERRQ(ierr);} 2055 ierr = PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, &qweights);CHKERRQ(ierr); 2056 /* For each fine grid cell */ 2057 for (cell = cStart; cell < cEnd; ++cell) { 2058 Vec pointVec; 2059 PetscScalar *pV; 2060 PetscSF coarseCellSF = NULL; 2061 const PetscSFNode *coarseCells; 2062 PetscInt numCoarseCells, cpdim, q, c, j; 2063 PetscInt *findices, *cindices; 2064 PetscInt numFIndices, numCIndices; 2065 2066 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2067 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 2068 /* Get points from the quadrature */ 2069 ierr = VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);CHKERRQ(ierr); 2070 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 2071 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2072 for (q = 0; q < Nq; ++q) { 2073 const PetscReal xi0[3] = {-1., -1., -1.}; 2074 2075 /* Transform point to real space */ 2076 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x); 2077 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 2078 } 2079 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2080 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 2081 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 2082 /* Update matrix */ 2083 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 2084 if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 2085 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2086 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 2087 PetscReal pVReal[3]; 2088 const PetscReal xi0[3] = {-1., -1., -1.}; 2089 2090 2091 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2092 /* Transform points from real space to coarse reference space */ 2093 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 2094 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 2095 CoordinatesRealToRef(dim, dim, xi0, v0c, invJc, pVReal, x); 2096 2097 if (id == PETSCFE_CLASSID) { 2098 PetscFE fe = (PetscFE) obj; 2099 PetscReal *B; 2100 2101 /* Evaluate coarse basis on contained point */ 2102 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 2103 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 2104 /* Get elemMat entries by multiplying by weight */ 2105 for (i = 0; i < numFIndices; ++i) { 2106 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 2107 for (j = 0; j < cpdim; ++j) { 2108 for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*Bfine[(ccell*numFIndices + i)*Nc + c]*qweights[ccell*Nc + c]*detJ; 2109 } 2110 /* Update interpolator */ 2111 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 2112 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 2113 ierr = MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);CHKERRQ(ierr); 2114 } 2115 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 2116 } else { 2117 cpdim = 1; 2118 for (i = 0; i < numFIndices; ++i) { 2119 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 2120 for (j = 0; j < cpdim; ++j) { 2121 for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*1.0*qweights[ccell*Nc + c]*detJ; 2122 } 2123 /* Update interpolator */ 2124 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 2125 ierr = PetscPrintf(PETSC_COMM_SELF, "Nq: %d %d Nf: %d %d Nc: %d %d\n", ccell, Nq, i, numFIndices, j, numCIndices);CHKERRQ(ierr); 2126 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 2127 ierr = MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);CHKERRQ(ierr); 2128 } 2129 } 2130 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2131 } 2132 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2133 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 2134 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 2135 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2136 } 2137 } 2138 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 2139 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 2140 ierr = PetscFree(elemMat);CHKERRQ(ierr); 2141 ierr = MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2142 ierr = MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2143 PetscFunctionReturn(0); 2144 } 2145 2146 /*@ 2147 DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns 2148 2149 Input Parameters: 2150 + dmc - The coarse mesh 2151 - dmf - The fine mesh 2152 - user - The user context 2153 2154 Output Parameter: 2155 . sc - The mapping 2156 2157 Level: developer 2158 2159 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 2160 @*/ 2161 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 2162 { 2163 PetscDS prob; 2164 PetscFE *feRef; 2165 PetscFV *fvRef; 2166 Vec fv, cv; 2167 IS fis, cis; 2168 PetscSection fsection, fglobalSection, csection, cglobalSection; 2169 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 2170 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m; 2171 PetscErrorCode ierr; 2172 2173 PetscFunctionBegin; 2174 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2175 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 2176 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 2177 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 2178 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 2179 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 2180 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 2181 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 2182 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 2183 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 2184 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 2185 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 2186 for (f = 0; f < Nf; ++f) { 2187 PetscObject obj; 2188 PetscClassId id; 2189 PetscInt fNb = 0, Nc = 0; 2190 2191 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 2192 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2193 if (id == PETSCFE_CLASSID) { 2194 PetscFE fe = (PetscFE) obj; 2195 2196 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 2197 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 2198 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 2199 } else if (id == PETSCFV_CLASSID) { 2200 PetscFV fv = (PetscFV) obj; 2201 PetscDualSpace Q; 2202 2203 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 2204 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 2205 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 2206 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 2207 } 2208 fTotDim += fNb*Nc; 2209 } 2210 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 2211 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 2212 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 2213 PetscFE feC; 2214 PetscFV fvC; 2215 PetscDualSpace QF, QC; 2216 PetscInt NcF, NcC, fpdim, cpdim; 2217 2218 if (feRef[field]) { 2219 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 2220 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 2221 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 2222 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 2223 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 2224 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 2225 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 2226 } else { 2227 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 2228 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 2229 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 2230 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 2231 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 2232 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 2233 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 2234 } 2235 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); 2236 for (c = 0; c < cpdim; ++c) { 2237 PetscQuadrature cfunc; 2238 const PetscReal *cqpoints; 2239 PetscInt NpC; 2240 PetscBool found = PETSC_FALSE; 2241 2242 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 2243 ierr = PetscQuadratureGetData(cfunc, NULL, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 2244 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 2245 for (f = 0; f < fpdim; ++f) { 2246 PetscQuadrature ffunc; 2247 const PetscReal *fqpoints; 2248 PetscReal sum = 0.0; 2249 PetscInt NpF, comp; 2250 2251 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 2252 ierr = PetscQuadratureGetData(ffunc, NULL, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 2253 if (NpC != NpF) continue; 2254 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 2255 if (sum > 1.0e-9) continue; 2256 for (comp = 0; comp < NcC; ++comp) { 2257 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 2258 } 2259 found = PETSC_TRUE; 2260 break; 2261 } 2262 if (!found) { 2263 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 2264 if (fvRef[field]) { 2265 PetscInt comp; 2266 for (comp = 0; comp < NcC; ++comp) { 2267 cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp; 2268 } 2269 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 2270 } 2271 } 2272 offsetC += cpdim*NcC; 2273 offsetF += fpdim*NcF; 2274 } 2275 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 2276 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 2277 2278 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 2279 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 2280 ierr = VecGetOwnershipRange(cv, &startC, &endC);CHKERRQ(ierr); 2281 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 2282 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 2283 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 2284 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 2285 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 2286 for (c = cStart; c < cEnd; ++c) { 2287 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 2288 for (d = 0; d < cTotDim; ++d) { 2289 if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue; 2290 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]]); 2291 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 2292 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 2293 } 2294 } 2295 ierr = PetscFree(cmap);CHKERRQ(ierr); 2296 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 2297 2298 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 2299 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 2300 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 2301 ierr = ISDestroy(&cis);CHKERRQ(ierr); 2302 ierr = ISDestroy(&fis);CHKERRQ(ierr); 2303 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 2304 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 2305 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2306 PetscFunctionReturn(0); 2307 } 2308