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