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