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