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, PETSC_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, PETSC_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*Nc, &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", 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(ierr); 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 /*@ 1360 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 1361 1362 Input Parameters: 1363 + dmf - The fine mesh 1364 . dmc - The coarse mesh 1365 - user - The user context 1366 1367 Output Parameter: 1368 . In - The interpolation matrix 1369 1370 Level: developer 1371 1372 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1373 @*/ 1374 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 1375 { 1376 DM_Plex *mesh = (DM_Plex *) dmf->data; 1377 const char *name = "Interpolator"; 1378 PetscDS prob; 1379 PetscSection fsection, csection, globalFSection, globalCSection; 1380 PetscHashJK ht; 1381 PetscLayout rLayout; 1382 PetscInt *dnz, *onz; 1383 PetscInt locRows, rStart, rEnd; 1384 PetscReal *x, *v0, *J, *invJ, detJ; 1385 PetscReal *v0c, *Jc, *invJc, detJc; 1386 PetscScalar *elemMat; 1387 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1388 PetscErrorCode ierr; 1389 1390 PetscFunctionBegin; 1391 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1392 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1393 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1394 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1395 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1396 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1397 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1398 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1399 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1400 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1401 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1402 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1403 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1404 ierr = PetscMalloc1(totDim, &elemMat);CHKERRQ(ierr); 1405 1406 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1407 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1408 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1409 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1410 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1411 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1412 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1413 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1414 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1415 for (field = 0; field < Nf; ++field) { 1416 PetscObject obj; 1417 PetscClassId id; 1418 PetscDualSpace Q = NULL; 1419 PetscQuadrature f; 1420 const PetscReal *qpoints; 1421 PetscInt Nc, Np, fpdim, i, d; 1422 1423 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1424 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1425 if (id == PETSCFE_CLASSID) { 1426 PetscFE fe = (PetscFE) obj; 1427 1428 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1429 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1430 } else if (id == PETSCFV_CLASSID) { 1431 PetscFV fv = (PetscFV) obj; 1432 1433 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1434 Nc = 1; 1435 } 1436 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1437 /* For each fine grid cell */ 1438 for (cell = cStart; cell < cEnd; ++cell) { 1439 PetscInt *findices, *cindices; 1440 PetscInt numFIndices, numCIndices; 1441 1442 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1443 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1444 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim); 1445 for (i = 0; i < fpdim; ++i) { 1446 Vec pointVec; 1447 PetscScalar *pV; 1448 PetscSF coarseCellSF = NULL; 1449 const PetscSFNode *coarseCells; 1450 PetscInt numCoarseCells, q, c; 1451 1452 /* Get points from the dual basis functional quadrature */ 1453 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1454 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1455 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1456 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1457 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1458 for (q = 0; q < Np; ++q) { 1459 /* Transform point to real space */ 1460 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1461 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1462 } 1463 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1464 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1465 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1466 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 1467 /* Update preallocation info */ 1468 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1469 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1470 { 1471 PetscHashJKKey key; 1472 PetscHashJKIter missing, iter; 1473 1474 key.j = findices[i]; 1475 if (key.j >= 0) { 1476 /* Get indices for coarse elements */ 1477 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1478 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1479 for (c = 0; c < numCIndices; ++c) { 1480 key.k = cindices[c]; 1481 if (key.k < 0) continue; 1482 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1483 if (missing) { 1484 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1485 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1486 else ++onz[key.j-rStart]; 1487 } 1488 } 1489 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1490 } 1491 } 1492 } 1493 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1494 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1495 } 1496 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1497 } 1498 } 1499 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1500 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1501 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1502 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1503 for (field = 0; field < Nf; ++field) { 1504 PetscObject obj; 1505 PetscClassId id; 1506 PetscDualSpace Q = NULL; 1507 PetscQuadrature f; 1508 const PetscReal *qpoints, *qweights; 1509 PetscInt Nc, qNc, Np, fpdim, i, d; 1510 1511 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1512 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1513 if (id == PETSCFE_CLASSID) { 1514 PetscFE fe = (PetscFE) obj; 1515 1516 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1517 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1518 } else if (id == PETSCFV_CLASSID) { 1519 PetscFV fv = (PetscFV) obj; 1520 1521 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1522 Nc = 1; 1523 } 1524 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1525 /* For each fine grid cell */ 1526 for (cell = cStart; cell < cEnd; ++cell) { 1527 PetscInt *findices, *cindices; 1528 PetscInt numFIndices, numCIndices; 1529 1530 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1531 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1532 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim); 1533 for (i = 0; i < fpdim; ++i) { 1534 Vec pointVec; 1535 PetscScalar *pV; 1536 PetscSF coarseCellSF = NULL; 1537 const PetscSFNode *coarseCells; 1538 PetscInt numCoarseCells, cpdim, q, c, j; 1539 1540 /* Get points from the dual basis functional quadrature */ 1541 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1542 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);CHKERRQ(ierr); 1543 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); 1544 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1545 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1546 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1547 for (q = 0; q < Np; ++q) { 1548 /* Transform point to real space */ 1549 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1550 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1551 } 1552 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1553 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1554 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1555 /* Update preallocation info */ 1556 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1557 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1558 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1559 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1560 PetscReal pVReal[3]; 1561 1562 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1563 /* Transform points from real space to coarse reference space */ 1564 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 1565 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 1566 CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x); 1567 1568 if (id == PETSCFE_CLASSID) { 1569 PetscFE fe = (PetscFE) obj; 1570 PetscReal *B; 1571 1572 /* Evaluate coarse basis on contained point */ 1573 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1574 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 1575 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 1576 /* Get elemMat entries by multiplying by weight */ 1577 for (j = 0; j < cpdim; ++j) { 1578 for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*qweights[ccell*qNc + c]; 1579 } 1580 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1581 } else { 1582 cpdim = 1; 1583 for (j = 0; j < cpdim; ++j) { 1584 for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c]; 1585 } 1586 } 1587 /* Update interpolator */ 1588 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 1589 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 1590 ierr = MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1591 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1592 } 1593 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1594 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1595 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1596 } 1597 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1598 } 1599 } 1600 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 1601 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 1602 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1603 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1604 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1605 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1606 PetscFunctionReturn(0); 1607 } 1608 1609 /*@ 1610 DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns 1611 1612 Input Parameters: 1613 + dmc - The coarse mesh 1614 - dmf - The fine mesh 1615 - user - The user context 1616 1617 Output Parameter: 1618 . sc - The mapping 1619 1620 Level: developer 1621 1622 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1623 @*/ 1624 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 1625 { 1626 PetscDS prob; 1627 PetscFE *feRef; 1628 PetscFV *fvRef; 1629 Vec fv, cv; 1630 IS fis, cis; 1631 PetscSection fsection, fglobalSection, csection, cglobalSection; 1632 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 1633 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m; 1634 PetscErrorCode ierr; 1635 1636 PetscFunctionBegin; 1637 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1638 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1639 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1640 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1641 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1642 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1643 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1644 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1645 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1646 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1647 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1648 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1649 for (f = 0; f < Nf; ++f) { 1650 PetscObject obj; 1651 PetscClassId id; 1652 PetscInt fNb = 0, Nc = 0; 1653 1654 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1655 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1656 if (id == PETSCFE_CLASSID) { 1657 PetscFE fe = (PetscFE) obj; 1658 1659 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1660 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 1661 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1662 } else if (id == PETSCFV_CLASSID) { 1663 PetscFV fv = (PetscFV) obj; 1664 PetscDualSpace Q; 1665 1666 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1667 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1668 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 1669 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1670 } 1671 fTotDim += fNb*Nc; 1672 } 1673 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1674 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 1675 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 1676 PetscFE feC; 1677 PetscFV fvC; 1678 PetscDualSpace QF, QC; 1679 PetscInt NcF, NcC, fpdim, cpdim; 1680 1681 if (feRef[field]) { 1682 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 1683 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 1684 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 1685 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 1686 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1687 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 1688 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1689 } else { 1690 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 1691 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 1692 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 1693 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 1694 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1695 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 1696 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1697 } 1698 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); 1699 for (c = 0; c < cpdim; ++c) { 1700 PetscQuadrature cfunc; 1701 const PetscReal *cqpoints; 1702 PetscInt NpC; 1703 PetscBool found = PETSC_FALSE; 1704 1705 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 1706 ierr = PetscQuadratureGetData(cfunc, NULL, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 1707 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 1708 for (f = 0; f < fpdim; ++f) { 1709 PetscQuadrature ffunc; 1710 const PetscReal *fqpoints; 1711 PetscReal sum = 0.0; 1712 PetscInt NpF, comp; 1713 1714 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 1715 ierr = PetscQuadratureGetData(ffunc, NULL, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 1716 if (NpC != NpF) continue; 1717 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 1718 if (sum > 1.0e-9) continue; 1719 for (comp = 0; comp < NcC; ++comp) { 1720 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 1721 } 1722 found = PETSC_TRUE; 1723 break; 1724 } 1725 if (!found) { 1726 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 1727 if (fvRef[field]) { 1728 PetscInt comp; 1729 for (comp = 0; comp < NcC; ++comp) { 1730 cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp; 1731 } 1732 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 1733 } 1734 } 1735 offsetC += cpdim*NcC; 1736 offsetF += fpdim*NcF; 1737 } 1738 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 1739 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1740 1741 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 1742 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 1743 ierr = VecGetOwnershipRange(cv, &startC, &endC);CHKERRQ(ierr); 1744 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 1745 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 1746 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 1747 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 1748 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 1749 for (c = cStart; c < cEnd; ++c) { 1750 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 1751 for (d = 0; d < cTotDim; ++d) { 1752 if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue; 1753 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]]); 1754 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 1755 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 1756 } 1757 } 1758 ierr = PetscFree(cmap);CHKERRQ(ierr); 1759 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 1760 1761 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 1762 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 1763 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 1764 ierr = ISDestroy(&cis);CHKERRQ(ierr); 1765 ierr = ISDestroy(&fis);CHKERRQ(ierr); 1766 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 1767 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 1768 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1769 PetscFunctionReturn(0); 1770 } 1771