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