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