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[], PetscInt, const PetscScalar[], 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[], PetscInt, const PetscScalar[], 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[], PetscInt, const PetscScalar[], 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 *constants, *lgrad; 988 PetscReal *lintegral; 989 PetscInt *uOff, *uOff_x, *aOff = NULL; 990 PetscInt dim, numConstants, 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 = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 1008 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1009 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1010 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1011 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1012 numCells = cEnd - cStart; 1013 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1014 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1015 if (dmAux) { 1016 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1017 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 1018 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1019 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1020 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 1021 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, NULL);CHKERRQ(ierr); 1022 } 1023 for (f = 0; f < Nf; ++f) { 1024 PetscObject obj; 1025 PetscClassId id; 1026 1027 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1028 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1029 if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;} 1030 } 1031 if (useFVM) { 1032 Vec grad; 1033 PetscInt fStart, fEnd; 1034 PetscBool compGrad; 1035 1036 ierr = PetscFVGetComputeGradients(fvm, &compGrad);CHKERRQ(ierr); 1037 ierr = PetscFVSetComputeGradients(fvm, PETSC_TRUE);CHKERRQ(ierr); 1038 ierr = DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);CHKERRQ(ierr); 1039 ierr = DMPlexComputeGradientFVM(dm, fvm, faceGeometryFVM, cellGeometryFVM, &dmGrad);CHKERRQ(ierr); 1040 ierr = PetscFVSetComputeGradients(fvm, compGrad);CHKERRQ(ierr); 1041 ierr = VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);CHKERRQ(ierr); 1042 ierr = VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1043 /* Reconstruct and limit cell gradients */ 1044 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1045 ierr = DMGetGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1046 ierr = DMPlexReconstructGradients_Internal(dm, fvm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, localX, grad);CHKERRQ(ierr); 1047 /* Communicate gradient values */ 1048 ierr = DMGetLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1049 ierr = DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1050 ierr = DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1051 ierr = DMRestoreGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1052 /* Handle non-essential (e.g. outflow) boundary values */ 1053 ierr = DMPlexInsertBoundaryValues(dm, PETSC_FALSE, localX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);CHKERRQ(ierr); 1054 ierr = VecGetArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1055 } 1056 ierr = PetscMalloc3(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeomFEM);CHKERRQ(ierr); 1057 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1058 for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;} 1059 for (c = cStart; c < cEnd; ++c) { 1060 PetscScalar *x = NULL; 1061 PetscInt i; 1062 1063 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeomFEM[c].v0, cgeomFEM[c].J, cgeomFEM[c].invJ, &cgeomFEM[c].detJ);CHKERRQ(ierr); 1064 if (cgeomFEM[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeomFEM[c].detJ, c); 1065 ierr = DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1066 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 1067 ierr = DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1068 if (dmAux) { 1069 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1070 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 1071 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1072 } 1073 } 1074 for (f = 0; f < Nf; ++f) { 1075 PetscObject obj; 1076 PetscClassId id; 1077 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 1078 1079 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1080 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1081 if (id == PETSCFE_CLASSID) { 1082 PetscFE fe = (PetscFE) obj; 1083 PetscQuadrature q; 1084 PetscInt Nq, Nb; 1085 1086 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1087 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1088 ierr = PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 1089 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1090 blockSize = Nb*Nq; 1091 batchSize = numBlocks * blockSize; 1092 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1093 numChunks = numCells / (numBatches*batchSize); 1094 Ne = numChunks*numBatches*batchSize; 1095 Nr = numCells % (numBatches*batchSize); 1096 offset = numCells - Nr; 1097 ierr = PetscFEIntegrate(fe, prob, f, Ne, cgeomFEM, u, probAux, a, lintegral);CHKERRQ(ierr); 1098 ierr = PetscFEIntegrate(fe, prob, f, Nr, &cgeomFEM[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);CHKERRQ(ierr); 1099 } else if (id == PETSCFV_CLASSID) { 1100 /* PetscFV fv = (PetscFV) obj; */ 1101 PetscInt foff; 1102 PetscPointFunc obj_func; 1103 PetscScalar lint; 1104 1105 ierr = PetscDSGetObjective(prob, f, &obj_func);CHKERRQ(ierr); 1106 ierr = PetscDSGetFieldOffset(prob, f, &foff);CHKERRQ(ierr); 1107 if (obj_func) { 1108 for (c = 0; c < numCells; ++c) { 1109 PetscScalar *u_x; 1110 1111 ierr = DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);CHKERRQ(ierr); 1112 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); 1113 lintegral[f] += PetscRealPart(lint)*cgeomFVM[c].volume; 1114 } 1115 } 1116 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 1117 } 1118 if (useFVM) { 1119 ierr = VecRestoreArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1120 ierr = VecRestoreArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);CHKERRQ(ierr); 1121 ierr = VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1122 ierr = DMRestoreLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1123 ierr = VecDestroy(&faceGeometryFVM);CHKERRQ(ierr); 1124 ierr = VecDestroy(&cellGeometryFVM);CHKERRQ(ierr); 1125 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr); 1126 } 1127 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1128 if (mesh->printFEM) { 1129 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");CHKERRQ(ierr); 1130 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);CHKERRQ(ierr);} 1131 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 1132 } 1133 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1134 ierr = MPIU_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr); 1135 ierr = PetscFree3(lintegral,u,cgeomFEM);CHKERRQ(ierr); 1136 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1137 PetscFunctionReturn(0); 1138 } 1139 1140 /*@ 1141 DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1142 1143 Input Parameters: 1144 + dmf - The fine mesh 1145 . dmc - The coarse mesh 1146 - user - The user context 1147 1148 Output Parameter: 1149 . In - The interpolation matrix 1150 1151 Level: developer 1152 1153 .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1154 @*/ 1155 PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user) 1156 { 1157 DM_Plex *mesh = (DM_Plex *) dmc->data; 1158 const char *name = "Interpolator"; 1159 PetscDS prob; 1160 PetscFE *feRef; 1161 PetscFV *fvRef; 1162 PetscSection fsection, fglobalSection; 1163 PetscSection csection, cglobalSection; 1164 PetscScalar *elemMat; 1165 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c; 1166 PetscInt cTotDim, rTotDim = 0; 1167 PetscErrorCode ierr; 1168 1169 PetscFunctionBegin; 1170 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1171 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1172 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1173 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1174 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1175 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1176 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1177 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1178 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1179 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1180 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr); 1181 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1182 for (f = 0; f < Nf; ++f) { 1183 PetscObject obj; 1184 PetscClassId id; 1185 PetscInt rNb = 0, Nc = 0; 1186 1187 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1188 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1189 if (id == PETSCFE_CLASSID) { 1190 PetscFE fe = (PetscFE) obj; 1191 1192 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1193 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1194 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1195 } else if (id == PETSCFV_CLASSID) { 1196 PetscFV fv = (PetscFV) obj; 1197 PetscDualSpace Q; 1198 1199 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1200 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1201 ierr = PetscDualSpaceGetDimension(Q, &rNb);CHKERRQ(ierr); 1202 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1203 } 1204 rTotDim += rNb; 1205 } 1206 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1207 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 1208 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 1209 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1210 PetscDualSpace Qref; 1211 PetscQuadrature f; 1212 const PetscReal *qpoints, *qweights; 1213 PetscReal *points; 1214 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1215 1216 /* Compose points from all dual basis functionals */ 1217 if (feRef[fieldI]) { 1218 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1219 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 1220 } else { 1221 ierr = PetscFVGetDualSpace(fvRef[fieldI], &Qref);CHKERRQ(ierr); 1222 ierr = PetscFVGetNumComponents(fvRef[fieldI], &Nc);CHKERRQ(ierr); 1223 } 1224 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1225 for (i = 0; i < fpdim; ++i) { 1226 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1227 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1228 npoints += Np; 1229 } 1230 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1231 for (i = 0, k = 0; i < fpdim; ++i) { 1232 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1233 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1234 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1235 } 1236 1237 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1238 PetscObject obj; 1239 PetscClassId id; 1240 PetscReal *B; 1241 PetscInt NcJ = 0, cpdim = 0, j, qNc; 1242 1243 ierr = PetscDSGetDiscretization(prob, fieldJ, &obj);CHKERRQ(ierr); 1244 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1245 if (id == PETSCFE_CLASSID) { 1246 PetscFE fe = (PetscFE) obj; 1247 1248 /* Evaluate basis at points */ 1249 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 1250 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1251 /* For now, fields only interpolate themselves */ 1252 if (fieldI == fieldJ) { 1253 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); 1254 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1255 for (i = 0, k = 0; i < fpdim; ++i) { 1256 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1257 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);CHKERRQ(ierr); 1258 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); 1259 for (p = 0; p < Np; ++p, ++k) { 1260 for (j = 0; j < cpdim; ++j) { 1261 /* NOTE: This is not quite right, unless fpdim == number of fine grid functional quad points */ 1262 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p*qNc+c]; 1263 } 1264 } 1265 } 1266 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1267 } 1268 } else if (id == PETSCFV_CLASSID) { 1269 PetscFV fv = (PetscFV) obj; 1270 1271 /* Evaluate constant function at points */ 1272 ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr); 1273 cpdim = 1; 1274 /* For now, fields only interpolate themselves */ 1275 if (fieldI == fieldJ) { 1276 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); 1277 for (i = 0, k = 0; i < fpdim; ++i) { 1278 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1279 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);CHKERRQ(ierr); 1280 if (qNc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in quadrature %D does not match coarse field %D", qNc, NcJ); 1281 for (p = 0; p < Np; ++p, ++k) { 1282 for (j = 0; j < cpdim; ++j) { 1283 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p*qNc+c]; 1284 } 1285 } 1286 } 1287 } 1288 } 1289 offsetJ += cpdim*NcJ; 1290 } 1291 offsetI += fpdim*Nc; 1292 ierr = PetscFree(points);CHKERRQ(ierr); 1293 } 1294 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1295 /* Preallocate matrix */ 1296 { 1297 Mat preallocator; 1298 PetscScalar *vals; 1299 PetscInt *cellCIndices, *cellFIndices; 1300 PetscInt locRows, locCols, cell; 1301 1302 ierr = MatGetLocalSize(In, &locRows, &locCols);CHKERRQ(ierr); 1303 ierr = MatCreate(PetscObjectComm((PetscObject) In), &preallocator);CHKERRQ(ierr); 1304 ierr = MatSetType(preallocator, MATPREALLOCATOR);CHKERRQ(ierr); 1305 ierr = MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1306 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 1307 ierr = PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1308 for (cell = cStart; cell < cEnd; ++cell) { 1309 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1310 ierr = MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);CHKERRQ(ierr); 1311 } 1312 ierr = PetscFree3(vals,cellCIndices,cellFIndices);CHKERRQ(ierr); 1313 ierr = MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1314 ierr = MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1315 ierr = MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);CHKERRQ(ierr); 1316 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 1317 } 1318 /* Fill matrix */ 1319 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1320 for (c = cStart; c < cEnd; ++c) { 1321 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1322 } 1323 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1324 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1325 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1326 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1327 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1328 if (mesh->printFEM) { 1329 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1330 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1331 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1332 } 1333 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1334 PetscFunctionReturn(0); 1335 } 1336 1337 /*@ 1338 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 1339 1340 Input Parameters: 1341 + dmf - The fine mesh 1342 . dmc - The coarse mesh 1343 - user - The user context 1344 1345 Output Parameter: 1346 . In - The interpolation matrix 1347 1348 Level: developer 1349 1350 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1351 @*/ 1352 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 1353 { 1354 DM_Plex *mesh = (DM_Plex *) dmf->data; 1355 const char *name = "Interpolator"; 1356 PetscDS prob; 1357 PetscSection fsection, csection, globalFSection, globalCSection; 1358 PetscHashJK ht; 1359 PetscLayout rLayout; 1360 PetscInt *dnz, *onz; 1361 PetscInt locRows, rStart, rEnd; 1362 PetscReal *x, *v0, *J, *invJ, detJ; 1363 PetscReal *v0c, *Jc, *invJc, detJc; 1364 PetscScalar *elemMat; 1365 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1366 PetscErrorCode ierr; 1367 1368 PetscFunctionBegin; 1369 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1370 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1371 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1372 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1373 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1374 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1375 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1376 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1377 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1378 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1379 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1380 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1381 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1382 ierr = PetscMalloc1(totDim, &elemMat);CHKERRQ(ierr); 1383 1384 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1385 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1386 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1387 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1388 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1389 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1390 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1391 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1392 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1393 for (field = 0; field < Nf; ++field) { 1394 PetscObject obj; 1395 PetscClassId id; 1396 PetscDualSpace Q = NULL; 1397 PetscQuadrature f; 1398 const PetscReal *qpoints; 1399 PetscInt Nc, Np, fpdim, i, d; 1400 1401 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1402 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1403 if (id == PETSCFE_CLASSID) { 1404 PetscFE fe = (PetscFE) obj; 1405 1406 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1407 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1408 } else if (id == PETSCFV_CLASSID) { 1409 PetscFV fv = (PetscFV) obj; 1410 1411 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1412 Nc = 1; 1413 } 1414 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1415 /* For each fine grid cell */ 1416 for (cell = cStart; cell < cEnd; ++cell) { 1417 PetscInt *findices, *cindices; 1418 PetscInt numFIndices, numCIndices; 1419 1420 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1421 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1422 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim); 1423 for (i = 0; i < fpdim; ++i) { 1424 Vec pointVec; 1425 PetscScalar *pV; 1426 PetscSF coarseCellSF = NULL; 1427 const PetscSFNode *coarseCells; 1428 PetscInt numCoarseCells, q, c; 1429 1430 /* Get points from the dual basis functional quadrature */ 1431 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1432 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1433 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1434 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1435 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1436 for (q = 0; q < Np; ++q) { 1437 /* Transform point to real space */ 1438 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1439 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1440 } 1441 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1442 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1443 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1444 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 1445 /* Update preallocation info */ 1446 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1447 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1448 { 1449 PetscHashJKKey key; 1450 PetscHashJKIter missing, iter; 1451 1452 key.j = findices[i]; 1453 if (key.j >= 0) { 1454 /* Get indices for coarse elements */ 1455 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1456 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1457 for (c = 0; c < numCIndices; ++c) { 1458 key.k = cindices[c]; 1459 if (key.k < 0) continue; 1460 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1461 if (missing) { 1462 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1463 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1464 else ++onz[key.j-rStart]; 1465 } 1466 } 1467 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1468 } 1469 } 1470 } 1471 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1472 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1473 } 1474 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1475 } 1476 } 1477 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1478 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1479 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1480 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1481 for (field = 0; field < Nf; ++field) { 1482 PetscObject obj; 1483 PetscClassId id; 1484 PetscDualSpace Q = NULL; 1485 PetscQuadrature f; 1486 const PetscReal *qpoints, *qweights; 1487 PetscInt Nc, qNc, Np, fpdim, i, d; 1488 1489 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1490 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1491 if (id == PETSCFE_CLASSID) { 1492 PetscFE fe = (PetscFE) obj; 1493 1494 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1495 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1496 } else if (id == PETSCFV_CLASSID) { 1497 PetscFV fv = (PetscFV) obj; 1498 1499 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1500 Nc = 1; 1501 } 1502 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1503 /* For each fine grid cell */ 1504 for (cell = cStart; cell < cEnd; ++cell) { 1505 PetscInt *findices, *cindices; 1506 PetscInt numFIndices, numCIndices; 1507 1508 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1509 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1510 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim); 1511 for (i = 0; i < fpdim; ++i) { 1512 Vec pointVec; 1513 PetscScalar *pV; 1514 PetscSF coarseCellSF = NULL; 1515 const PetscSFNode *coarseCells; 1516 PetscInt numCoarseCells, cpdim, q, c, j; 1517 1518 /* Get points from the dual basis functional quadrature */ 1519 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1520 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);CHKERRQ(ierr); 1521 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); 1522 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1523 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1524 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1525 for (q = 0; q < Np; ++q) { 1526 /* Transform point to real space */ 1527 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1528 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1529 } 1530 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1531 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1532 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1533 /* Update preallocation info */ 1534 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1535 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1536 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1537 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1538 PetscReal pVReal[3]; 1539 1540 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1541 /* Transform points from real space to coarse reference space */ 1542 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 1543 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 1544 CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x); 1545 1546 if (id == PETSCFE_CLASSID) { 1547 PetscFE fe = (PetscFE) obj; 1548 PetscReal *B; 1549 1550 /* Evaluate coarse basis on contained point */ 1551 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1552 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 1553 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 1554 /* Get elemMat entries by multiplying by weight */ 1555 for (j = 0; j < cpdim; ++j) { 1556 for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*qweights[ccell*qNc + c]; 1557 } 1558 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1559 } else { 1560 cpdim = 1; 1561 for (j = 0; j < cpdim; ++j) { 1562 for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c]; 1563 } 1564 } 1565 /* Update interpolator */ 1566 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 1567 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 1568 ierr = MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1569 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1570 } 1571 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1572 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1573 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1574 } 1575 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1576 } 1577 } 1578 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 1579 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 1580 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1581 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1582 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1583 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1584 PetscFunctionReturn(0); 1585 } 1586 1587 /*@ 1588 DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns 1589 1590 Input Parameters: 1591 + dmc - The coarse mesh 1592 - dmf - The fine mesh 1593 - user - The user context 1594 1595 Output Parameter: 1596 . sc - The mapping 1597 1598 Level: developer 1599 1600 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1601 @*/ 1602 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 1603 { 1604 PetscDS prob; 1605 PetscFE *feRef; 1606 PetscFV *fvRef; 1607 Vec fv, cv; 1608 IS fis, cis; 1609 PetscSection fsection, fglobalSection, csection, cglobalSection; 1610 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 1611 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m; 1612 PetscErrorCode ierr; 1613 1614 PetscFunctionBegin; 1615 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1616 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1617 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1618 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1619 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1620 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1621 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1622 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1623 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1624 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1625 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1626 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1627 for (f = 0; f < Nf; ++f) { 1628 PetscObject obj; 1629 PetscClassId id; 1630 PetscInt fNb = 0, Nc = 0; 1631 1632 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1633 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1634 if (id == PETSCFE_CLASSID) { 1635 PetscFE fe = (PetscFE) obj; 1636 1637 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1638 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 1639 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1640 } else if (id == PETSCFV_CLASSID) { 1641 PetscFV fv = (PetscFV) obj; 1642 PetscDualSpace Q; 1643 1644 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1645 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1646 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 1647 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1648 } 1649 fTotDim += fNb*Nc; 1650 } 1651 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1652 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 1653 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 1654 PetscFE feC; 1655 PetscFV fvC; 1656 PetscDualSpace QF, QC; 1657 PetscInt NcF, NcC, fpdim, cpdim; 1658 1659 if (feRef[field]) { 1660 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 1661 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 1662 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 1663 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 1664 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1665 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 1666 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1667 } else { 1668 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 1669 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 1670 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 1671 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 1672 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1673 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 1674 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1675 } 1676 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); 1677 for (c = 0; c < cpdim; ++c) { 1678 PetscQuadrature cfunc; 1679 const PetscReal *cqpoints; 1680 PetscInt NpC; 1681 PetscBool found = PETSC_FALSE; 1682 1683 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 1684 ierr = PetscQuadratureGetData(cfunc, NULL, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 1685 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 1686 for (f = 0; f < fpdim; ++f) { 1687 PetscQuadrature ffunc; 1688 const PetscReal *fqpoints; 1689 PetscReal sum = 0.0; 1690 PetscInt NpF, comp; 1691 1692 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 1693 ierr = PetscQuadratureGetData(ffunc, NULL, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 1694 if (NpC != NpF) continue; 1695 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 1696 if (sum > 1.0e-9) continue; 1697 for (comp = 0; comp < NcC; ++comp) { 1698 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 1699 } 1700 found = PETSC_TRUE; 1701 break; 1702 } 1703 if (!found) { 1704 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 1705 if (fvRef[field]) { 1706 PetscInt comp; 1707 for (comp = 0; comp < NcC; ++comp) { 1708 cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp; 1709 } 1710 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 1711 } 1712 } 1713 offsetC += cpdim*NcC; 1714 offsetF += fpdim*NcF; 1715 } 1716 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 1717 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1718 1719 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 1720 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 1721 ierr = VecGetOwnershipRange(cv, &startC, &endC);CHKERRQ(ierr); 1722 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 1723 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 1724 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 1725 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 1726 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 1727 for (c = cStart; c < cEnd; ++c) { 1728 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 1729 for (d = 0; d < cTotDim; ++d) { 1730 if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue; 1731 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]]); 1732 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 1733 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 1734 } 1735 } 1736 ierr = PetscFree(cmap);CHKERRQ(ierr); 1737 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 1738 1739 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 1740 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 1741 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 1742 ierr = ISDestroy(&cis);CHKERRQ(ierr); 1743 ierr = ISDestroy(&fis);CHKERRQ(ierr); 1744 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 1745 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 1746 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1747 PetscFunctionReturn(0); 1748 } 1749