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