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