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 const PetscInt debug = 0; 569 PetscSection section; 570 PetscQuadrature quad; 571 Vec localX; 572 PetscScalar *funcVal, *interpolant; 573 PetscReal *coords, *detJ, *J; 574 PetscReal localDiff = 0.0; 575 const PetscReal *quadWeights; 576 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset; 577 PetscErrorCode ierr; 578 579 PetscFunctionBegin; 580 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 581 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 582 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 583 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 584 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 585 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 586 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 587 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 588 for (field = 0; field < numFields; ++field) { 589 PetscObject obj; 590 PetscClassId id; 591 PetscInt Nc; 592 593 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 594 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 595 if (id == PETSCFE_CLASSID) { 596 PetscFE fe = (PetscFE) obj; 597 598 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 599 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 600 } else if (id == PETSCFV_CLASSID) { 601 PetscFV fv = (PetscFV) obj; 602 603 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 604 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 605 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 606 numComponents += Nc; 607 } 608 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);CHKERRQ(ierr); 609 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 610 ierr = PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);CHKERRQ(ierr); 611 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 612 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 613 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 614 for (c = cStart; c < cEnd; ++c) { 615 PetscScalar *x = NULL; 616 PetscReal elemDiff = 0.0; 617 PetscInt qc = 0; 618 619 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);CHKERRQ(ierr); 620 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 621 622 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 623 PetscObject obj; 624 PetscClassId id; 625 void * const ctx = ctxs ? ctxs[field] : NULL; 626 PetscInt Nb, Nc, q, fc; 627 628 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 629 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 630 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 631 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 632 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 633 if (debug) { 634 char title[1024]; 635 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 636 ierr = DMPrintCellVector(c, title, Nb, &x[fieldOffset]);CHKERRQ(ierr); 637 } 638 for (q = 0; q < Nq; ++q) { 639 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); 640 ierr = (*funcs[field])(coordDim, time, &coords[coordDim * q], Nc, funcVal, ctx); 641 if (ierr) { 642 PetscErrorCode ierr2; 643 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 644 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 645 ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2); 646 CHKERRQ(ierr); 647 } 648 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 649 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 650 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 651 for (fc = 0; fc < Nc; ++fc) { 652 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 653 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);} 654 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]; 655 } 656 } 657 fieldOffset += Nb; 658 qc += Nc; 659 } 660 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 661 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 662 localDiff += elemDiff; 663 } 664 ierr = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr); 665 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 666 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 667 *diff = PetscSqrtReal(*diff); 668 PetscFunctionReturn(0); 669 } 670 671 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) 672 { 673 const PetscInt debug = 0; 674 PetscSection section; 675 PetscQuadrature quad; 676 Vec localX; 677 PetscScalar *funcVal, *interpolant; 678 const PetscReal *quadPoints, *quadWeights; 679 PetscReal *coords, *realSpaceDer, *J, *invJ, *detJ; 680 PetscReal localDiff = 0.0; 681 PetscInt dim, coordDim, qNc = 0, Nq = 0, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset; 682 PetscErrorCode ierr; 683 684 PetscFunctionBegin; 685 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 686 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 687 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 688 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 689 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 690 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 691 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 692 for (field = 0; field < numFields; ++field) { 693 PetscFE fe; 694 PetscInt Nc; 695 696 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 697 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 698 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 699 numComponents += Nc; 700 } 701 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 702 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 703 /* ierr = DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 704 ierr = PetscMalloc7(numComponents,&funcVal,coordDim*Nq,&coords,coordDim*Nq,&realSpaceDer,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ,numComponents,&interpolant,Nq,&detJ);CHKERRQ(ierr); 705 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 706 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 707 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 708 for (c = cStart; c < cEnd; ++c) { 709 PetscScalar *x = NULL; 710 PetscReal elemDiff = 0.0; 711 PetscInt qc = 0; 712 713 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, invJ, detJ);CHKERRQ(ierr); 714 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 715 716 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 717 PetscFE fe; 718 void * const ctx = ctxs ? ctxs[field] : NULL; 719 PetscReal *basisDer; 720 PetscInt Nb, Nc, q, fc; 721 722 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 723 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 724 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 725 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 726 if (debug) { 727 char title[1024]; 728 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 729 ierr = DMPrintCellVector(c, title, Nb, &x[fieldOffset]);CHKERRQ(ierr); 730 } 731 for (q = 0; q < Nq; ++q) { 732 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); 733 ierr = (*funcs[field])(coordDim, time, &coords[q*coordDim], n, numFields, funcVal, ctx); 734 if (ierr) { 735 PetscErrorCode ierr2; 736 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 737 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 738 ierr2 = PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);CHKERRQ(ierr2); 739 CHKERRQ(ierr); 740 } 741 ierr = PetscFEInterpolateGradient_Static(fe, &x[fieldOffset], coordDim, invJ, n, q, interpolant);CHKERRQ(ierr); 742 for (fc = 0; fc < Nc; ++fc) { 743 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 744 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);} 745 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]; 746 } 747 } 748 fieldOffset += Nb; 749 qc += Nc; 750 } 751 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 752 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 753 localDiff += elemDiff; 754 } 755 ierr = PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);CHKERRQ(ierr); 756 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 757 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 758 *diff = PetscSqrtReal(*diff); 759 PetscFunctionReturn(0); 760 } 761 762 PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 763 { 764 const PetscInt debug = 0; 765 PetscSection section; 766 PetscQuadrature quad; 767 Vec localX; 768 PetscScalar *funcVal, *interpolant; 769 PetscReal *coords, *detJ, *J; 770 PetscReal *localDiff; 771 const PetscReal *quadPoints, *quadWeights; 772 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset; 773 PetscErrorCode ierr; 774 775 PetscFunctionBegin; 776 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 777 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 778 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 779 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 780 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 781 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 782 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 783 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 784 for (field = 0; field < numFields; ++field) { 785 PetscObject obj; 786 PetscClassId id; 787 PetscInt Nc; 788 789 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 790 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 791 if (id == PETSCFE_CLASSID) { 792 PetscFE fe = (PetscFE) obj; 793 794 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 795 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 796 } else if (id == PETSCFV_CLASSID) { 797 PetscFV fv = (PetscFV) obj; 798 799 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 800 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 801 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 802 numComponents += Nc; 803 } 804 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 805 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 806 ierr = PetscCalloc6(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);CHKERRQ(ierr); 807 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 808 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 809 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 810 for (c = cStart; c < cEnd; ++c) { 811 PetscScalar *x = NULL; 812 PetscInt qc = 0; 813 814 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);CHKERRQ(ierr); 815 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 816 817 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 818 PetscObject obj; 819 PetscClassId id; 820 void * const ctx = ctxs ? ctxs[field] : NULL; 821 PetscInt Nb, Nc, q, fc; 822 823 PetscReal elemDiff = 0.0; 824 825 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 826 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 827 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 828 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 829 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 830 if (debug) { 831 char title[1024]; 832 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 833 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 834 } 835 for (q = 0; q < Nq; ++q) { 836 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); 837 ierr = (*funcs[field])(coordDim, time, &coords[coordDim*q], numFields, funcVal, ctx); 838 if (ierr) { 839 PetscErrorCode ierr2; 840 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 841 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 842 ierr2 = PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2); 843 CHKERRQ(ierr); 844 } 845 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 846 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 847 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 848 for (fc = 0; fc < Nc; ++fc) { 849 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 850 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);} 851 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]; 852 } 853 } 854 fieldOffset += Nb; 855 qc += Nc; 856 localDiff[field] += elemDiff; 857 } 858 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 859 } 860 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 861 ierr = MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 862 for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]); 863 ierr = PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr); 864 PetscFunctionReturn(0); 865 } 866 867 /*@C 868 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. 869 870 Input Parameters: 871 + dm - The DM 872 . time - The time 873 . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation 874 . ctxs - Optional array of contexts to pass to each function, or NULL. 875 - X - The coefficient vector u_h 876 877 Output Parameter: 878 . D - A Vec which holds the difference ||u - u_h||_2 for each cell 879 880 Level: developer 881 882 .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff() 883 @*/ 884 PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D) 885 { 886 PetscSection section; 887 PetscQuadrature quad; 888 Vec localX; 889 PetscScalar *funcVal, *interpolant; 890 PetscReal *coords, *detJ, *J; 891 const PetscReal *quadPoints, *quadWeights; 892 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset; 893 PetscErrorCode ierr; 894 895 PetscFunctionBegin; 896 ierr = VecSet(D, 0.0);CHKERRQ(ierr); 897 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 898 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 899 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 900 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 901 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 902 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 903 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 904 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 905 for (field = 0; field < numFields; ++field) { 906 PetscObject obj; 907 PetscClassId id; 908 PetscInt Nc; 909 910 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 911 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 912 if (id == PETSCFE_CLASSID) { 913 PetscFE fe = (PetscFE) obj; 914 915 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 916 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 917 } else if (id == PETSCFV_CLASSID) { 918 PetscFV fv = (PetscFV) obj; 919 920 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 921 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 922 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 923 numComponents += Nc; 924 } 925 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 926 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 927 ierr = PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);CHKERRQ(ierr); 928 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 929 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 930 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 931 for (c = cStart; c < cEnd; ++c) { 932 PetscScalar *x = NULL; 933 PetscScalar elemDiff = 0.0; 934 PetscInt qc = 0; 935 936 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);CHKERRQ(ierr); 937 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 938 939 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 940 PetscObject obj; 941 PetscClassId id; 942 void * const ctx = ctxs ? ctxs[field] : NULL; 943 PetscInt Nb, Nc, q, fc; 944 945 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 946 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 947 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 948 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 949 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 950 if (funcs[field]) { 951 for (q = 0; q < Nq; ++q) { 952 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); 953 ierr = (*funcs[field])(dim, time, &coords[q*coordDim], Nc, funcVal, ctx); 954 if (ierr) { 955 PetscErrorCode ierr2; 956 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 957 ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2); 958 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 959 CHKERRQ(ierr); 960 } 961 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 962 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 963 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 964 for (fc = 0; fc < Nc; ++fc) { 965 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 966 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]; 967 } 968 } 969 } 970 fieldOffset += Nb; 971 qc += Nc; 972 } 973 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 974 ierr = VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);CHKERRQ(ierr); 975 } 976 ierr = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr); 977 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 978 ierr = VecSqrtAbs(D);CHKERRQ(ierr); 979 PetscFunctionReturn(0); 980 } 981 982 /*@ 983 DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user 984 985 Input Parameters: 986 + dm - The mesh 987 . X - Local input vector 988 - user - The user context 989 990 Output Parameter: 991 . integral - Local integral for each field 992 993 Level: developer 994 995 .seealso: DMPlexComputeResidualFEM() 996 @*/ 997 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user) 998 { 999 DM_Plex *mesh = (DM_Plex *) dm->data; 1000 DM dmAux, dmGrad; 1001 Vec localX, A, cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL; 1002 PetscDS prob, probAux = NULL; 1003 PetscSection section, sectionAux; 1004 PetscFV fvm = NULL; 1005 PetscFECellGeom *cgeomFEM; 1006 PetscFVFaceGeom *fgeomFVM; 1007 PetscFVCellGeom *cgeomFVM; 1008 PetscScalar *u, *a = NULL; 1009 const PetscScalar *constants, *lgrad; 1010 PetscReal *lintegral; 1011 PetscInt *uOff, *uOff_x, *aOff = NULL; 1012 PetscInt dim, numConstants, Nf, NfAux = 0, f, numCells, cStart, cEnd, cEndInterior, c; 1013 PetscInt totDim, totDimAux; 1014 PetscBool useFVM = PETSC_FALSE; 1015 PetscErrorCode ierr; 1016 1017 PetscFunctionBegin; 1018 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1019 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1020 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 1021 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1022 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1023 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1024 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1025 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1026 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1027 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 1028 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 1029 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 1030 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1031 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1032 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1033 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1034 numCells = cEnd - cStart; 1035 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1036 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1037 if (dmAux) { 1038 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1039 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 1040 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1041 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1042 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 1043 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, NULL);CHKERRQ(ierr); 1044 } 1045 for (f = 0; f < Nf; ++f) { 1046 PetscObject obj; 1047 PetscClassId id; 1048 1049 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1050 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1051 if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;} 1052 } 1053 if (useFVM) { 1054 Vec grad; 1055 PetscInt fStart, fEnd; 1056 PetscBool compGrad; 1057 1058 ierr = PetscFVGetComputeGradients(fvm, &compGrad);CHKERRQ(ierr); 1059 ierr = PetscFVSetComputeGradients(fvm, PETSC_TRUE);CHKERRQ(ierr); 1060 ierr = DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);CHKERRQ(ierr); 1061 ierr = DMPlexComputeGradientFVM(dm, fvm, faceGeometryFVM, cellGeometryFVM, &dmGrad);CHKERRQ(ierr); 1062 ierr = PetscFVSetComputeGradients(fvm, compGrad);CHKERRQ(ierr); 1063 ierr = VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);CHKERRQ(ierr); 1064 ierr = VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1065 /* Reconstruct and limit cell gradients */ 1066 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1067 ierr = DMGetGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1068 ierr = DMPlexReconstructGradients_Internal(dm, fvm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, localX, grad);CHKERRQ(ierr); 1069 /* Communicate gradient values */ 1070 ierr = DMGetLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1071 ierr = DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1072 ierr = DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1073 ierr = DMRestoreGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1074 /* Handle non-essential (e.g. outflow) boundary values */ 1075 ierr = DMPlexInsertBoundaryValues(dm, PETSC_FALSE, localX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);CHKERRQ(ierr); 1076 ierr = VecGetArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1077 } 1078 ierr = PetscMalloc3(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeomFEM);CHKERRQ(ierr); 1079 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1080 for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;} 1081 for (c = cStart; c < cEnd; ++c) { 1082 PetscScalar *x = NULL; 1083 PetscInt i; 1084 1085 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeomFEM[c].v0, cgeomFEM[c].J, cgeomFEM[c].invJ, &cgeomFEM[c].detJ);CHKERRQ(ierr); 1086 if (cgeomFEM[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeomFEM[c].detJ, c); 1087 ierr = DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1088 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 1089 ierr = DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1090 if (dmAux) { 1091 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1092 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 1093 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1094 } 1095 } 1096 for (f = 0; f < Nf; ++f) { 1097 PetscObject obj; 1098 PetscClassId id; 1099 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 1100 1101 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1102 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1103 if (id == PETSCFE_CLASSID) { 1104 PetscFE fe = (PetscFE) obj; 1105 PetscQuadrature q; 1106 PetscInt Nq, Nb; 1107 1108 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1109 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1110 ierr = PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 1111 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1112 blockSize = Nb*Nq; 1113 batchSize = numBlocks * blockSize; 1114 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1115 numChunks = numCells / (numBatches*batchSize); 1116 Ne = numChunks*numBatches*batchSize; 1117 Nr = numCells % (numBatches*batchSize); 1118 offset = numCells - Nr; 1119 ierr = PetscFEIntegrate(fe, prob, f, Ne, cgeomFEM, u, probAux, a, lintegral);CHKERRQ(ierr); 1120 ierr = PetscFEIntegrate(fe, prob, f, Nr, &cgeomFEM[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);CHKERRQ(ierr); 1121 } else if (id == PETSCFV_CLASSID) { 1122 /* PetscFV fv = (PetscFV) obj; */ 1123 PetscInt foff; 1124 PetscPointFunc obj_func; 1125 PetscScalar lint; 1126 1127 ierr = PetscDSGetObjective(prob, f, &obj_func);CHKERRQ(ierr); 1128 ierr = PetscDSGetFieldOffset(prob, f, &foff);CHKERRQ(ierr); 1129 if (obj_func) { 1130 for (c = 0; c < numCells; ++c) { 1131 PetscScalar *u_x; 1132 1133 ierr = DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);CHKERRQ(ierr); 1134 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); 1135 lintegral[f] += PetscRealPart(lint)*cgeomFVM[c].volume; 1136 } 1137 } 1138 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 1139 } 1140 if (useFVM) { 1141 ierr = VecRestoreArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1142 ierr = VecRestoreArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);CHKERRQ(ierr); 1143 ierr = VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1144 ierr = DMRestoreLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1145 ierr = VecDestroy(&faceGeometryFVM);CHKERRQ(ierr); 1146 ierr = VecDestroy(&cellGeometryFVM);CHKERRQ(ierr); 1147 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr); 1148 } 1149 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1150 if (mesh->printFEM) { 1151 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");CHKERRQ(ierr); 1152 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);CHKERRQ(ierr);} 1153 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 1154 } 1155 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1156 ierr = MPIU_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr); 1157 ierr = PetscFree3(lintegral,u,cgeomFEM);CHKERRQ(ierr); 1158 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1159 PetscFunctionReturn(0); 1160 } 1161 1162 /*@ 1163 DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1164 1165 Input Parameters: 1166 + dmf - The fine mesh 1167 . dmc - The coarse mesh 1168 - user - The user context 1169 1170 Output Parameter: 1171 . In - The interpolation matrix 1172 1173 Level: developer 1174 1175 .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1176 @*/ 1177 PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user) 1178 { 1179 DM_Plex *mesh = (DM_Plex *) dmc->data; 1180 const char *name = "Interpolator"; 1181 PetscDS prob; 1182 PetscFE *feRef; 1183 PetscFV *fvRef; 1184 PetscSection fsection, fglobalSection; 1185 PetscSection csection, cglobalSection; 1186 PetscScalar *elemMat; 1187 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c; 1188 PetscInt cTotDim, rTotDim = 0; 1189 PetscErrorCode ierr; 1190 1191 PetscFunctionBegin; 1192 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1193 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1194 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1195 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1196 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1197 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1198 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1199 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1200 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1201 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1202 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr); 1203 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1204 for (f = 0; f < Nf; ++f) { 1205 PetscObject obj; 1206 PetscClassId id; 1207 PetscInt rNb = 0, Nc = 0; 1208 1209 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1210 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1211 if (id == PETSCFE_CLASSID) { 1212 PetscFE fe = (PetscFE) obj; 1213 1214 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1215 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1216 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1217 } else if (id == PETSCFV_CLASSID) { 1218 PetscFV fv = (PetscFV) obj; 1219 PetscDualSpace Q; 1220 1221 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1222 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1223 ierr = PetscDualSpaceGetDimension(Q, &rNb);CHKERRQ(ierr); 1224 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1225 } 1226 rTotDim += rNb; 1227 } 1228 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1229 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 1230 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 1231 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1232 PetscDualSpace Qref; 1233 PetscQuadrature f; 1234 const PetscReal *qpoints, *qweights; 1235 PetscReal *points; 1236 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1237 1238 /* Compose points from all dual basis functionals */ 1239 if (feRef[fieldI]) { 1240 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1241 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 1242 } else { 1243 ierr = PetscFVGetDualSpace(fvRef[fieldI], &Qref);CHKERRQ(ierr); 1244 ierr = PetscFVGetNumComponents(fvRef[fieldI], &Nc);CHKERRQ(ierr); 1245 } 1246 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1247 for (i = 0; i < fpdim; ++i) { 1248 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1249 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1250 npoints += Np; 1251 } 1252 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1253 for (i = 0, k = 0; i < fpdim; ++i) { 1254 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1255 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1256 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1257 } 1258 1259 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1260 PetscObject obj; 1261 PetscClassId id; 1262 PetscReal *B; 1263 PetscInt NcJ = 0, cpdim = 0, j, qNc; 1264 1265 ierr = PetscDSGetDiscretization(prob, fieldJ, &obj);CHKERRQ(ierr); 1266 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1267 if (id == PETSCFE_CLASSID) { 1268 PetscFE fe = (PetscFE) obj; 1269 1270 /* Evaluate basis at points */ 1271 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 1272 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1273 /* For now, fields only interpolate themselves */ 1274 if (fieldI == fieldJ) { 1275 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); 1276 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1277 for (i = 0, k = 0; i < fpdim; ++i) { 1278 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1279 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);CHKERRQ(ierr); 1280 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); 1281 for (p = 0; p < Np; ++p, ++k) { 1282 for (j = 0; j < cpdim; ++j) { 1283 /* NOTE: This is not quite right, unless fpdim == number of fine grid functional quad points */ 1284 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p*qNc+c]; 1285 } 1286 } 1287 } 1288 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1289 } 1290 } else if (id == PETSCFV_CLASSID) { 1291 PetscFV fv = (PetscFV) obj; 1292 1293 /* Evaluate constant function at points */ 1294 ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr); 1295 cpdim = 1; 1296 /* For now, fields only interpolate themselves */ 1297 if (fieldI == fieldJ) { 1298 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); 1299 for (i = 0, k = 0; i < fpdim; ++i) { 1300 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1301 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);CHKERRQ(ierr); 1302 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); 1303 for (p = 0; p < Np; ++p, ++k) { 1304 for (j = 0; j < cpdim; ++j) { 1305 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p*qNc+c]; 1306 } 1307 } 1308 } 1309 } 1310 } 1311 offsetJ += cpdim*NcJ; 1312 } 1313 offsetI += fpdim*Nc; 1314 ierr = PetscFree(points);CHKERRQ(ierr); 1315 } 1316 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1317 /* Preallocate matrix */ 1318 { 1319 Mat preallocator; 1320 PetscScalar *vals; 1321 PetscInt *cellCIndices, *cellFIndices; 1322 PetscInt locRows, locCols, cell; 1323 1324 ierr = MatGetLocalSize(In, &locRows, &locCols);CHKERRQ(ierr); 1325 ierr = MatCreate(PetscObjectComm((PetscObject) In), &preallocator);CHKERRQ(ierr); 1326 ierr = MatSetType(preallocator, MATPREALLOCATOR);CHKERRQ(ierr); 1327 ierr = MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1328 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 1329 ierr = PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1330 for (cell = cStart; cell < cEnd; ++cell) { 1331 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1332 ierr = MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);CHKERRQ(ierr); 1333 } 1334 ierr = PetscFree3(vals,cellCIndices,cellFIndices);CHKERRQ(ierr); 1335 ierr = MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1336 ierr = MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1337 ierr = MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);CHKERRQ(ierr); 1338 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 1339 } 1340 /* Fill matrix */ 1341 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1342 for (c = cStart; c < cEnd; ++c) { 1343 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1344 } 1345 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1346 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1347 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1348 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1349 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1350 if (mesh->printFEM) { 1351 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1352 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1353 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1354 } 1355 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1356 PetscFunctionReturn(0); 1357 } 1358 1359 PetscErrorCode DMPlexComputeMassMatrixNested(DM dmc, DM dmf, Mat mass, void *user) 1360 { 1361 SETERRQ(PetscObjectComm((PetscObject) dmc), PETSC_ERR_SUP, "Laziness"); 1362 } 1363 1364 /*@ 1365 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 1366 1367 Input Parameters: 1368 + dmf - The fine mesh 1369 . dmc - The coarse mesh 1370 - user - The user context 1371 1372 Output Parameter: 1373 . In - The interpolation matrix 1374 1375 Level: developer 1376 1377 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1378 @*/ 1379 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 1380 { 1381 DM_Plex *mesh = (DM_Plex *) dmf->data; 1382 const char *name = "Interpolator"; 1383 PetscDS prob; 1384 PetscSection fsection, csection, globalFSection, globalCSection; 1385 PetscHashJK ht; 1386 PetscLayout rLayout; 1387 PetscInt *dnz, *onz; 1388 PetscInt locRows, rStart, rEnd; 1389 PetscReal *x, *v0, *J, *invJ, detJ; 1390 PetscReal *v0c, *Jc, *invJc, detJc; 1391 PetscScalar *elemMat; 1392 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1393 PetscErrorCode ierr; 1394 1395 PetscFunctionBegin; 1396 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1397 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1398 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1399 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1400 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1401 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1402 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1403 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1404 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1405 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1406 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1407 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1408 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1409 ierr = PetscMalloc1(totDim, &elemMat);CHKERRQ(ierr); 1410 1411 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1412 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1413 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1414 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1415 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1416 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1417 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1418 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1419 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1420 for (field = 0; field < Nf; ++field) { 1421 PetscObject obj; 1422 PetscClassId id; 1423 PetscDualSpace Q = NULL; 1424 PetscQuadrature f; 1425 const PetscReal *qpoints; 1426 PetscInt Nc, Np, fpdim, i, d; 1427 1428 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1429 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1430 if (id == PETSCFE_CLASSID) { 1431 PetscFE fe = (PetscFE) obj; 1432 1433 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1434 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1435 } else if (id == PETSCFV_CLASSID) { 1436 PetscFV fv = (PetscFV) obj; 1437 1438 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1439 Nc = 1; 1440 } 1441 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1442 /* For each fine grid cell */ 1443 for (cell = cStart; cell < cEnd; ++cell) { 1444 PetscInt *findices, *cindices; 1445 PetscInt numFIndices, numCIndices; 1446 1447 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1448 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1449 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim); 1450 for (i = 0; i < fpdim; ++i) { 1451 Vec pointVec; 1452 PetscScalar *pV; 1453 PetscSF coarseCellSF = NULL; 1454 const PetscSFNode *coarseCells; 1455 PetscInt numCoarseCells, q, c; 1456 1457 /* Get points from the dual basis functional quadrature */ 1458 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1459 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1460 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1461 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1462 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1463 for (q = 0; q < Np; ++q) { 1464 /* Transform point to real space */ 1465 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1466 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1467 } 1468 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1469 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1470 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1471 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 1472 /* Update preallocation info */ 1473 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1474 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1475 { 1476 PetscHashJKKey key; 1477 PetscHashJKIter missing, iter; 1478 1479 key.j = findices[i]; 1480 if (key.j >= 0) { 1481 /* Get indices for coarse elements */ 1482 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1483 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1484 for (c = 0; c < numCIndices; ++c) { 1485 key.k = cindices[c]; 1486 if (key.k < 0) continue; 1487 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1488 if (missing) { 1489 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1490 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1491 else ++onz[key.j-rStart]; 1492 } 1493 } 1494 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1495 } 1496 } 1497 } 1498 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1499 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1500 } 1501 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1502 } 1503 } 1504 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1505 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1506 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1507 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1508 for (field = 0; field < Nf; ++field) { 1509 PetscObject obj; 1510 PetscClassId id; 1511 PetscDualSpace Q = NULL; 1512 PetscQuadrature f; 1513 const PetscReal *qpoints, *qweights; 1514 PetscInt Nc, qNc, Np, fpdim, i, d; 1515 1516 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1517 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1518 if (id == PETSCFE_CLASSID) { 1519 PetscFE fe = (PetscFE) obj; 1520 1521 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1522 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1523 } else if (id == PETSCFV_CLASSID) { 1524 PetscFV fv = (PetscFV) obj; 1525 1526 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1527 Nc = 1; 1528 } 1529 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1530 /* For each fine grid cell */ 1531 for (cell = cStart; cell < cEnd; ++cell) { 1532 PetscInt *findices, *cindices; 1533 PetscInt numFIndices, numCIndices; 1534 1535 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1536 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1537 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim); 1538 for (i = 0; i < fpdim; ++i) { 1539 Vec pointVec; 1540 PetscScalar *pV; 1541 PetscSF coarseCellSF = NULL; 1542 const PetscSFNode *coarseCells; 1543 PetscInt numCoarseCells, cpdim, q, c, j; 1544 1545 /* Get points from the dual basis functional quadrature */ 1546 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1547 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);CHKERRQ(ierr); 1548 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); 1549 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1550 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1551 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1552 for (q = 0; q < Np; ++q) { 1553 /* Transform point to real space */ 1554 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1555 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1556 } 1557 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1558 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1559 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1560 /* Update preallocation info */ 1561 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1562 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1563 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1564 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1565 PetscReal pVReal[3]; 1566 1567 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1568 /* Transform points from real space to coarse reference space */ 1569 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 1570 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 1571 CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x); 1572 1573 if (id == PETSCFE_CLASSID) { 1574 PetscFE fe = (PetscFE) obj; 1575 PetscReal *B; 1576 1577 /* Evaluate coarse basis on contained point */ 1578 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1579 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 1580 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 1581 /* Get elemMat entries by multiplying by weight */ 1582 for (j = 0; j < cpdim; ++j) { 1583 for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*qweights[ccell*qNc + c]; 1584 } 1585 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1586 } else { 1587 cpdim = 1; 1588 for (j = 0; j < cpdim; ++j) { 1589 for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c]; 1590 } 1591 } 1592 /* Update interpolator */ 1593 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 1594 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 1595 ierr = MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1596 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1597 } 1598 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1599 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1600 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1601 } 1602 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1603 } 1604 } 1605 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 1606 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 1607 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1608 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1609 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1610 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1611 PetscFunctionReturn(0); 1612 } 1613 1614 /*@ 1615 DMPlexComputeMassMatrixGeneral - Form the local portion of the mass matrix M from the coarse DM to a non-nested fine DM. 1616 1617 Input Parameters: 1618 + dmf - The fine mesh 1619 . dmc - The coarse mesh 1620 - user - The user context 1621 1622 Output Parameter: 1623 . mass - The mass matrix 1624 1625 Level: developer 1626 1627 .seealso: DMPlexComputeMassMatrixNested(), DMPlexComputeInterpolatorNested(), DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1628 @*/ 1629 PetscErrorCode DMPlexComputeMassMatrixGeneral(DM dmc, DM dmf, Mat mass, void *user) 1630 { 1631 DM_Plex *mesh = (DM_Plex *) dmf->data; 1632 const char *name = "Mass Matrix"; 1633 PetscDS prob; 1634 PetscSection fsection, csection, globalFSection, globalCSection; 1635 PetscHashJK ht; 1636 PetscLayout rLayout; 1637 PetscInt *dnz, *onz; 1638 PetscInt locRows, rStart, rEnd; 1639 PetscReal *x, *v0, *J, *invJ, detJ; 1640 PetscReal *v0c, *Jc, *invJc, detJc; 1641 PetscScalar *elemMat; 1642 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1643 PetscErrorCode ierr; 1644 1645 PetscFunctionBegin; 1646 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1647 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1648 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1649 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1650 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1651 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1652 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1653 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1654 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1655 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1656 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1657 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1658 ierr = PetscMalloc1(totDim, &elemMat);CHKERRQ(ierr); 1659 1660 ierr = MatGetLocalSize(mass, &locRows, NULL);CHKERRQ(ierr); 1661 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) mass), &rLayout);CHKERRQ(ierr); 1662 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1663 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1664 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1665 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1666 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1667 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1668 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1669 for (field = 0; field < Nf; ++field) { 1670 PetscObject obj; 1671 PetscClassId id; 1672 PetscQuadrature quad; 1673 const PetscReal *qpoints; 1674 PetscInt Nq, Nc, i, d; 1675 1676 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1677 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1678 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetQuadrature((PetscFE) obj, &quad);CHKERRQ(ierr);} 1679 else {ierr = PetscFVGetQuadrature((PetscFV) obj, &quad);CHKERRQ(ierr);} 1680 ierr = PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, NULL);CHKERRQ(ierr); 1681 /* For each fine grid cell */ 1682 for (cell = cStart; cell < cEnd; ++cell) { 1683 Vec pointVec; 1684 PetscScalar *pV; 1685 PetscSF coarseCellSF = NULL; 1686 const PetscSFNode *coarseCells; 1687 PetscInt numCoarseCells, q, c; 1688 PetscInt *findices, *cindices; 1689 PetscInt numFIndices, numCIndices; 1690 1691 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1692 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1693 /* Get points from the quadrature */ 1694 ierr = VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);CHKERRQ(ierr); 1695 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1696 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1697 for (q = 0; q < Nq; ++q) { 1698 /* Transform point to real space */ 1699 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1700 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1701 } 1702 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1703 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1704 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1705 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 1706 /* Update preallocation info */ 1707 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1708 if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1709 { 1710 PetscHashJKKey key; 1711 PetscHashJKIter missing, iter; 1712 1713 for (i = 0; i < numFIndices; ++i) { 1714 key.j = findices[i]; 1715 if (key.j >= 0) { 1716 /* Get indices for coarse elements */ 1717 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1718 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1719 for (c = 0; c < numCIndices; ++c) { 1720 key.k = cindices[c]; 1721 if (key.k < 0) continue; 1722 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1723 if (missing) { 1724 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1725 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1726 else ++onz[key.j-rStart]; 1727 } 1728 } 1729 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1730 } 1731 } 1732 } 1733 } 1734 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1735 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1736 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1737 } 1738 } 1739 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1740 ierr = MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1741 ierr = MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1742 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1743 for (field = 0; field < Nf; ++field) { 1744 PetscObject obj; 1745 PetscClassId id; 1746 PetscQuadrature quad; 1747 PetscReal *Bfine; 1748 const PetscReal *qpoints, *qweights; 1749 PetscInt Nq, Nc, i, d; 1750 1751 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1752 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1753 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetQuadrature((PetscFE) obj, &quad);CHKERRQ(ierr);ierr = PetscFEGetDefaultTabulation((PetscFE) obj, &Bfine, NULL, NULL);CHKERRQ(ierr);} 1754 else {ierr = PetscFVGetQuadrature((PetscFV) obj, &quad);CHKERRQ(ierr);} 1755 ierr = PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, &qweights);CHKERRQ(ierr); 1756 /* For each fine grid cell */ 1757 for (cell = cStart; cell < cEnd; ++cell) { 1758 Vec pointVec; 1759 PetscScalar *pV; 1760 PetscSF coarseCellSF = NULL; 1761 const PetscSFNode *coarseCells; 1762 PetscInt numCoarseCells, cpdim, q, c, j; 1763 PetscInt *findices, *cindices; 1764 PetscInt numFIndices, numCIndices; 1765 1766 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1767 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1768 /* Get points from the quadrature */ 1769 ierr = VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);CHKERRQ(ierr); 1770 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1771 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1772 for (q = 0; q < Nq; ++q) { 1773 /* Transform point to real space */ 1774 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1775 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1776 } 1777 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1778 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1779 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1780 /* Update matrix */ 1781 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1782 if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1783 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1784 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1785 PetscReal pVReal[3]; 1786 1787 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1788 /* Transform points from real space to coarse reference space */ 1789 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 1790 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 1791 CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x); 1792 1793 if (id == PETSCFE_CLASSID) { 1794 PetscFE fe = (PetscFE) obj; 1795 PetscReal *B; 1796 1797 /* Evaluate coarse basis on contained point */ 1798 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1799 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 1800 /* Get elemMat entries by multiplying by weight */ 1801 for (i = 0; i < numFIndices; ++i) { 1802 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 1803 for (j = 0; j < cpdim; ++j) { 1804 for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*Bfine[(ccell*numFIndices + i)*Nc + c]*qweights[ccell*Nc + c]*detJ; 1805 } 1806 /* Update interpolator */ 1807 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 1808 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 1809 ierr = MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);CHKERRQ(ierr); 1810 } 1811 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1812 } else { 1813 cpdim = 1; 1814 for (i = 0; i < numFIndices; ++i) { 1815 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 1816 for (j = 0; j < cpdim; ++j) { 1817 for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*1.0*qweights[ccell*Nc + c]*detJ; 1818 } 1819 /* Update interpolator */ 1820 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 1821 ierr = PetscPrintf(PETSC_COMM_SELF, "Nq: %d %d Nf: %d %d Nc: %d %d\n", ccell, Nq, i, numFIndices, j, numCIndices);CHKERRQ(ierr); 1822 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 1823 ierr = MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);CHKERRQ(ierr); 1824 } 1825 } 1826 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1827 } 1828 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1829 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1830 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1831 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1832 } 1833 } 1834 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 1835 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 1836 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1837 ierr = MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1838 ierr = MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1839 PetscFunctionReturn(0); 1840 } 1841 1842 /*@ 1843 DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns 1844 1845 Input Parameters: 1846 + dmc - The coarse mesh 1847 - dmf - The fine mesh 1848 - user - The user context 1849 1850 Output Parameter: 1851 . sc - The mapping 1852 1853 Level: developer 1854 1855 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1856 @*/ 1857 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 1858 { 1859 PetscDS prob; 1860 PetscFE *feRef; 1861 PetscFV *fvRef; 1862 Vec fv, cv; 1863 IS fis, cis; 1864 PetscSection fsection, fglobalSection, csection, cglobalSection; 1865 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 1866 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m; 1867 PetscErrorCode ierr; 1868 1869 PetscFunctionBegin; 1870 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1871 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1872 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1873 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1874 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1875 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1876 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1877 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1878 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1879 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1880 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1881 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1882 for (f = 0; f < Nf; ++f) { 1883 PetscObject obj; 1884 PetscClassId id; 1885 PetscInt fNb = 0, Nc = 0; 1886 1887 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1888 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1889 if (id == PETSCFE_CLASSID) { 1890 PetscFE fe = (PetscFE) obj; 1891 1892 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1893 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 1894 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1895 } else if (id == PETSCFV_CLASSID) { 1896 PetscFV fv = (PetscFV) obj; 1897 PetscDualSpace Q; 1898 1899 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1900 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1901 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 1902 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1903 } 1904 fTotDim += fNb*Nc; 1905 } 1906 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1907 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 1908 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 1909 PetscFE feC; 1910 PetscFV fvC; 1911 PetscDualSpace QF, QC; 1912 PetscInt NcF, NcC, fpdim, cpdim; 1913 1914 if (feRef[field]) { 1915 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 1916 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 1917 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 1918 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 1919 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1920 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 1921 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1922 } else { 1923 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 1924 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 1925 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 1926 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 1927 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1928 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 1929 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1930 } 1931 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); 1932 for (c = 0; c < cpdim; ++c) { 1933 PetscQuadrature cfunc; 1934 const PetscReal *cqpoints; 1935 PetscInt NpC; 1936 PetscBool found = PETSC_FALSE; 1937 1938 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 1939 ierr = PetscQuadratureGetData(cfunc, NULL, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 1940 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 1941 for (f = 0; f < fpdim; ++f) { 1942 PetscQuadrature ffunc; 1943 const PetscReal *fqpoints; 1944 PetscReal sum = 0.0; 1945 PetscInt NpF, comp; 1946 1947 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 1948 ierr = PetscQuadratureGetData(ffunc, NULL, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 1949 if (NpC != NpF) continue; 1950 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 1951 if (sum > 1.0e-9) continue; 1952 for (comp = 0; comp < NcC; ++comp) { 1953 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 1954 } 1955 found = PETSC_TRUE; 1956 break; 1957 } 1958 if (!found) { 1959 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 1960 if (fvRef[field]) { 1961 PetscInt comp; 1962 for (comp = 0; comp < NcC; ++comp) { 1963 cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp; 1964 } 1965 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 1966 } 1967 } 1968 offsetC += cpdim*NcC; 1969 offsetF += fpdim*NcF; 1970 } 1971 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 1972 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1973 1974 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 1975 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 1976 ierr = VecGetOwnershipRange(cv, &startC, &endC);CHKERRQ(ierr); 1977 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 1978 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 1979 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 1980 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 1981 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 1982 for (c = cStart; c < cEnd; ++c) { 1983 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 1984 for (d = 0; d < cTotDim; ++d) { 1985 if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue; 1986 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]]); 1987 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 1988 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 1989 } 1990 } 1991 ierr = PetscFree(cmap);CHKERRQ(ierr); 1992 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 1993 1994 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 1995 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 1996 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 1997 ierr = ISDestroy(&cis);CHKERRQ(ierr); 1998 ierr = ISDestroy(&fis);CHKERRQ(ierr); 1999 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 2000 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 2001 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2002 PetscFunctionReturn(0); 2003 } 2004