1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petscsf.h> 3 4 #include <petsc/private/petscfeimpl.h> 5 #include <petsc/private/petscfvimpl.h> 6 7 /*@ 8 DMPlexGetScale - Get the scale for the specified fundamental unit 9 10 Not collective 11 12 Input Arguments: 13 + dm - the DM 14 - unit - The SI unit 15 16 Output Argument: 17 . scale - The value used to scale all quantities with this unit 18 19 Level: advanced 20 21 .seealso: DMPlexSetScale(), PetscUnit 22 @*/ 23 PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale) 24 { 25 DM_Plex *mesh = (DM_Plex*) dm->data; 26 27 PetscFunctionBegin; 28 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 29 PetscValidPointer(scale, 3); 30 *scale = mesh->scale[unit]; 31 PetscFunctionReturn(0); 32 } 33 34 /*@ 35 DMPlexSetScale - Set the scale for the specified fundamental unit 36 37 Not collective 38 39 Input Arguments: 40 + dm - the DM 41 . unit - The SI unit 42 - scale - The value used to scale all quantities with this unit 43 44 Level: advanced 45 46 .seealso: DMPlexGetScale(), PetscUnit 47 @*/ 48 PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale) 49 { 50 DM_Plex *mesh = (DM_Plex*) dm->data; 51 52 PetscFunctionBegin; 53 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 54 mesh->scale[unit] = scale; 55 PetscFunctionReturn(0); 56 } 57 58 PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k) 59 { 60 switch (i) { 61 case 0: 62 switch (j) { 63 case 0: return 0; 64 case 1: 65 switch (k) { 66 case 0: return 0; 67 case 1: return 0; 68 case 2: return 1; 69 } 70 case 2: 71 switch (k) { 72 case 0: return 0; 73 case 1: return -1; 74 case 2: return 0; 75 } 76 } 77 case 1: 78 switch (j) { 79 case 0: 80 switch (k) { 81 case 0: return 0; 82 case 1: return 0; 83 case 2: return -1; 84 } 85 case 1: return 0; 86 case 2: 87 switch (k) { 88 case 0: return 1; 89 case 1: return 0; 90 case 2: return 0; 91 } 92 } 93 case 2: 94 switch (j) { 95 case 0: 96 switch (k) { 97 case 0: return 0; 98 case 1: return 1; 99 case 2: return 0; 100 } 101 case 1: 102 switch (k) { 103 case 0: return -1; 104 case 1: return 0; 105 case 2: return 0; 106 } 107 case 2: return 0; 108 } 109 } 110 return 0; 111 } 112 113 static PetscErrorCode DMPlexProjectRigidBody_Private(PetscInt dim, PetscReal t, const PetscReal X[], PetscInt Nf, PetscScalar *mode, void *ctx) 114 { 115 PetscInt *ctxInt = (PetscInt *) ctx; 116 PetscInt dim2 = ctxInt[0]; 117 PetscInt d = ctxInt[1]; 118 PetscInt i, j, k = dim > 2 ? d - dim : d; 119 120 PetscFunctionBegin; 121 if (dim != dim2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input dimension %d does not match context dimension %d", dim, dim2); 122 for (i = 0; i < dim; i++) mode[i] = 0.; 123 if (d < dim) { 124 mode[d] = 1.; 125 } else { 126 for (i = 0; i < dim; i++) { 127 for (j = 0; j < dim; j++) { 128 mode[j] += epsilon(i, j, k)*X[i]; 129 } 130 } 131 } 132 PetscFunctionReturn(0); 133 } 134 135 /*@C 136 DMPlexCreateRigidBody - for the default global section, create rigid body modes from coordinates 137 138 Collective on DM 139 140 Input Arguments: 141 . dm - the DM 142 143 Output Argument: 144 . sp - the null space 145 146 Note: This is necessary to take account of Dirichlet conditions on the displacements 147 148 Level: advanced 149 150 .seealso: MatNullSpaceCreate() 151 @*/ 152 PetscErrorCode DMPlexCreateRigidBody(DM dm, MatNullSpace *sp) 153 { 154 MPI_Comm comm; 155 Vec mode[6]; 156 PetscSection section, globalSection; 157 PetscInt dim, dimEmbed, n, m, d, i, j; 158 PetscErrorCode ierr; 159 160 PetscFunctionBegin; 161 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 162 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 163 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 164 if (dim == 1) { 165 ierr = MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);CHKERRQ(ierr); 166 PetscFunctionReturn(0); 167 } 168 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 169 ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr); 170 ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr); 171 m = (dim*(dim+1))/2; 172 ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr); 173 ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr); 174 ierr = VecSetUp(mode[0]);CHKERRQ(ierr); 175 for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);} 176 for (d = 0; d < m; d++) { 177 PetscInt ctx[2]; 178 PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private; 179 void *voidctx = (void *) (&ctx[0]); 180 181 ctx[0] = dimEmbed; 182 ctx[1] = d; 183 ierr = DMProjectFunction(dm, 0.0, &func, &voidctx, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 184 } 185 for (i = 0; i < dim; ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);} 186 /* Orthonormalize system */ 187 for (i = dim; i < m; ++i) { 188 PetscScalar dots[6]; 189 190 ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr); 191 for (j = 0; j < i; ++j) dots[j] *= -1.0; 192 ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr); 193 ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr); 194 } 195 ierr = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);CHKERRQ(ierr); 196 for (i = 0; i< m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);} 197 PetscFunctionReturn(0); 198 } 199 200 /*@ 201 DMPlexSetMaxProjectionHeight - In DMPlexProjectXXXLocal() functions, the projected values of a basis function's dofs 202 are computed by associating the basis function with one of the mesh points in its transitively-closed support, and 203 evaluating the dual space basis of that point. A basis function is associated with the point in its 204 transitively-closed support whose mesh height is highest (w.r.t. DAG height), but not greater than the maximum 205 projection height, which is set with this function. By default, the maximum projection height is zero, which means 206 that only mesh cells are used to project basis functions. A height of one, for example, evaluates a cell-interior 207 basis functions using its cells dual space basis, but all other basis functions with the dual space basis of a face. 208 209 Input Parameters: 210 + dm - the DMPlex object 211 - height - the maximum projection height >= 0 212 213 Level: advanced 214 215 .seealso: DMPlexGetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal() 216 @*/ 217 PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height) 218 { 219 DM_Plex *plex = (DM_Plex *) dm->data; 220 221 PetscFunctionBegin; 222 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 223 plex->maxProjectionHeight = height; 224 PetscFunctionReturn(0); 225 } 226 227 /*@ 228 DMPlexGetMaxProjectionHeight - Get the maximum height (w.r.t. DAG) of mesh points used to evaluate dual bases in 229 DMPlexProjectXXXLocal() functions. 230 231 Input Parameters: 232 . dm - the DMPlex object 233 234 Output Parameters: 235 . height - the maximum projection height 236 237 Level: intermediate 238 239 .seealso: DMPlexSetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal() 240 @*/ 241 PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height) 242 { 243 DM_Plex *plex = (DM_Plex *) dm->data; 244 245 PetscFunctionBegin; 246 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 247 *height = plex->maxProjectionHeight; 248 PetscFunctionReturn(0); 249 } 250 251 /*@C 252 DMPlexInsertBoundaryValuesEssential - Insert boundary values into a local vector 253 254 Input Parameters: 255 + dm - The DM, with a PetscDS that matches the problem being constrained 256 . time - The time 257 . field - The field to constrain 258 . label - The DMLabel defining constrained points 259 . numids - The number of DMLabel ids for constrained points 260 . ids - An array of ids for constrained points 261 . func - A pointwise function giving boundary values 262 - ctx - An optional user context for bcFunc 263 264 Output Parameter: 265 . locX - A local vector to receives the boundary values 266 267 Level: developer 268 269 .seealso: DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary() 270 @*/ 271 PetscErrorCode DMPlexInsertBoundaryValuesEssential(DM dm, PetscReal time, PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void *ctx, Vec locX) 272 { 273 PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx); 274 void **ctxs; 275 PetscInt numFields; 276 PetscErrorCode ierr; 277 278 PetscFunctionBegin; 279 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 280 ierr = PetscCalloc2(numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 281 funcs[field] = func; 282 ctxs[field] = ctx; 283 ierr = DMProjectFunctionLabelLocal(dm, time, label, numids, ids, funcs, ctxs, INSERT_BC_VALUES, locX);CHKERRQ(ierr); 284 ierr = PetscFree2(funcs,ctxs);CHKERRQ(ierr); 285 PetscFunctionReturn(0); 286 } 287 288 /*@C 289 DMPlexInsertBoundaryValuesEssentialField - Insert boundary values into a local vector 290 291 Input Parameters: 292 + dm - The DM, with a PetscDS that matches the problem being constrained 293 . time - The time 294 . locU - A local vector with the input solution values 295 . field - The field to constrain 296 . label - The DMLabel defining constrained points 297 . numids - The number of DMLabel ids for constrained points 298 . ids - An array of ids for constrained points 299 . func - A pointwise function giving boundary values 300 - ctx - An optional user context for bcFunc 301 302 Output Parameter: 303 . locX - A local vector to receives the boundary values 304 305 Level: developer 306 307 .seealso: DMPlexInsertBoundaryValuesEssential(), DMAddBoundary() 308 @*/ 309 PetscErrorCode DMPlexInsertBoundaryValuesEssentialField(DM dm, PetscReal time, Vec locU, PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], 310 void (*func)(PetscInt, PetscInt, PetscInt, 311 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 312 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 313 PetscReal, const PetscReal[], 314 PetscScalar[]), 315 void *ctx, Vec locX) 316 { 317 void (**funcs)(PetscInt, PetscInt, PetscInt, 318 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 319 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 320 PetscReal, const PetscReal[], PetscScalar[]); 321 void **ctxs; 322 PetscInt numFields; 323 PetscErrorCode ierr; 324 325 PetscFunctionBegin; 326 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 327 ierr = PetscCalloc2(numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 328 funcs[field] = func; 329 ctxs[field] = ctx; 330 ierr = DMProjectFieldLabelLocal(dm, time, label, numids, ids, locU, funcs, INSERT_BC_VALUES, locX);CHKERRQ(ierr); 331 ierr = PetscFree2(funcs,ctxs);CHKERRQ(ierr); 332 PetscFunctionReturn(0); 333 } 334 335 /*@C 336 DMPlexInsertBoundaryValuesRiemann - Insert boundary values into a local vector 337 338 Input Parameters: 339 + dm - The DM, with a PetscDS that matches the problem being constrained 340 . time - The time 341 . faceGeometry - A vector with the FVM face geometry information 342 . cellGeometry - A vector with the FVM cell geometry information 343 . Grad - A vector with the FVM cell gradient information 344 . field - The field to constrain 345 . label - The DMLabel defining constrained points 346 . numids - The number of DMLabel ids for constrained points 347 . ids - An array of ids for constrained points 348 . func - A pointwise function giving boundary values 349 - ctx - An optional user context for bcFunc 350 351 Output Parameter: 352 . locX - A local vector to receives the boundary values 353 354 Note: This implementation currently ignores the numcomps/comps argument from DMAddBoundary() 355 356 Level: developer 357 358 .seealso: DMPlexInsertBoundaryValuesEssential(), DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary() 359 @*/ 360 PetscErrorCode DMPlexInsertBoundaryValuesRiemann(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad, PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], 361 PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX) 362 { 363 PetscDS prob; 364 PetscSF sf; 365 DM dmFace, dmCell, dmGrad; 366 const PetscScalar *facegeom, *cellgeom = NULL, *grad; 367 const PetscInt *leaves; 368 PetscScalar *x, *fx; 369 PetscInt dim, nleaves, loc, fStart, fEnd, pdim, i; 370 PetscErrorCode ierr, ierru = 0; 371 372 PetscFunctionBegin; 373 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 374 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 375 nleaves = PetscMax(0, nleaves); 376 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 377 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 378 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 379 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 380 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 381 if (cellGeometry) { 382 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 383 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 384 } 385 if (Grad) { 386 PetscFV fv; 387 388 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);CHKERRQ(ierr); 389 ierr = VecGetDM(Grad, &dmGrad);CHKERRQ(ierr); 390 ierr = VecGetArrayRead(Grad, &grad);CHKERRQ(ierr); 391 ierr = PetscFVGetNumComponents(fv, &pdim);CHKERRQ(ierr); 392 ierr = DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 393 } 394 ierr = VecGetArray(locX, &x);CHKERRQ(ierr); 395 for (i = 0; i < numids; ++i) { 396 IS faceIS; 397 const PetscInt *faces; 398 PetscInt numFaces, f; 399 400 ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr); 401 if (!faceIS) continue; /* No points with that id on this process */ 402 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 403 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 404 for (f = 0; f < numFaces; ++f) { 405 const PetscInt face = faces[f], *cells; 406 PetscFVFaceGeom *fg; 407 408 if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */ 409 ierr = PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);CHKERRQ(ierr); 410 if (loc >= 0) continue; 411 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 412 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 413 if (Grad) { 414 PetscFVCellGeom *cg; 415 PetscScalar *cx, *cgrad; 416 PetscScalar *xG; 417 PetscReal dx[3]; 418 PetscInt d; 419 420 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);CHKERRQ(ierr); 421 ierr = DMPlexPointLocalRead(dm, cells[0], x, &cx);CHKERRQ(ierr); 422 ierr = DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);CHKERRQ(ierr); 423 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 424 DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx); 425 for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx); 426 ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx); 427 if (ierru) { 428 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 429 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 430 goto cleanup; 431 } 432 } else { 433 PetscScalar *xI; 434 PetscScalar *xG; 435 436 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr); 437 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 438 ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx); 439 if (ierru) { 440 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 441 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 442 goto cleanup; 443 } 444 } 445 } 446 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 447 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 448 } 449 cleanup: 450 ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr); 451 if (Grad) { 452 ierr = DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 453 ierr = VecRestoreArrayRead(Grad, &grad);CHKERRQ(ierr); 454 } 455 if (cellGeometry) {ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr);} 456 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 457 CHKERRQ(ierru); 458 PetscFunctionReturn(0); 459 } 460 461 /*@ 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[], 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 *lgrad; 979 PetscReal *lintegral; 980 PetscInt *uOff, *uOff_x, *aOff = NULL; 981 PetscInt dim, 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 = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 999 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1000 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1001 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1002 numCells = cEnd - cStart; 1003 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1004 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1005 if (dmAux) { 1006 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1007 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 1008 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1009 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1010 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 1011 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, NULL);CHKERRQ(ierr); 1012 } 1013 for (f = 0; f < Nf; ++f) { 1014 PetscObject obj; 1015 PetscClassId id; 1016 1017 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1018 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1019 if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;} 1020 } 1021 if (useFVM) { 1022 Vec grad; 1023 PetscInt fStart, fEnd; 1024 PetscBool compGrad; 1025 1026 ierr = PetscFVGetComputeGradients(fvm, &compGrad);CHKERRQ(ierr); 1027 ierr = PetscFVSetComputeGradients(fvm, PETSC_TRUE);CHKERRQ(ierr); 1028 ierr = DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);CHKERRQ(ierr); 1029 ierr = DMPlexComputeGradientFVM(dm, fvm, faceGeometryFVM, cellGeometryFVM, &dmGrad);CHKERRQ(ierr); 1030 ierr = PetscFVSetComputeGradients(fvm, compGrad);CHKERRQ(ierr); 1031 ierr = VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);CHKERRQ(ierr); 1032 ierr = VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1033 /* Reconstruct and limit cell gradients */ 1034 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1035 ierr = DMGetGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1036 ierr = DMPlexReconstructGradients_Internal(dm, fvm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, localX, grad);CHKERRQ(ierr); 1037 /* Communicate gradient values */ 1038 ierr = DMGetLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1039 ierr = DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1040 ierr = DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1041 ierr = DMRestoreGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1042 /* Handle non-essential (e.g. outflow) boundary values */ 1043 ierr = DMPlexInsertBoundaryValues(dm, PETSC_FALSE, localX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);CHKERRQ(ierr); 1044 ierr = VecGetArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1045 } 1046 ierr = PetscMalloc3(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeomFEM);CHKERRQ(ierr); 1047 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1048 for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;} 1049 for (c = cStart; c < cEnd; ++c) { 1050 PetscScalar *x = NULL; 1051 PetscInt i; 1052 1053 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeomFEM[c].v0, cgeomFEM[c].J, cgeomFEM[c].invJ, &cgeomFEM[c].detJ);CHKERRQ(ierr); 1054 if (cgeomFEM[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeomFEM[c].detJ, c); 1055 ierr = DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1056 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 1057 ierr = DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1058 if (dmAux) { 1059 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1060 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 1061 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1062 } 1063 } 1064 for (f = 0; f < Nf; ++f) { 1065 PetscObject obj; 1066 PetscClassId id; 1067 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 1068 1069 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1070 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1071 if (id == PETSCFE_CLASSID) { 1072 PetscFE fe = (PetscFE) obj; 1073 PetscQuadrature q; 1074 PetscInt Nq, Nb; 1075 1076 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1077 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1078 ierr = PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 1079 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1080 blockSize = Nb*Nq; 1081 batchSize = numBlocks * blockSize; 1082 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1083 numChunks = numCells / (numBatches*batchSize); 1084 Ne = numChunks*numBatches*batchSize; 1085 Nr = numCells % (numBatches*batchSize); 1086 offset = numCells - Nr; 1087 ierr = PetscFEIntegrate(fe, prob, f, Ne, cgeomFEM, u, probAux, a, lintegral);CHKERRQ(ierr); 1088 ierr = PetscFEIntegrate(fe, prob, f, Nr, &cgeomFEM[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);CHKERRQ(ierr); 1089 } else if (id == PETSCFV_CLASSID) { 1090 /* PetscFV fv = (PetscFV) obj; */ 1091 PetscInt foff; 1092 PetscPointFunc obj_func; 1093 PetscScalar lint; 1094 1095 ierr = PetscDSGetObjective(prob, f, &obj_func);CHKERRQ(ierr); 1096 ierr = PetscDSGetFieldOffset(prob, f, &foff);CHKERRQ(ierr); 1097 if (obj_func) { 1098 for (c = 0; c < numCells; ++c) { 1099 PetscScalar *u_x; 1100 1101 ierr = DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);CHKERRQ(ierr); 1102 obj_func(dim, Nf, NfAux, uOff, uOff_x, &u[totDim*c+foff], NULL, u_x, aOff, NULL, a, NULL, NULL, 0.0, cgeomFVM[c].centroid, &lint); 1103 lintegral[f] += PetscRealPart(lint)*cgeomFVM[c].volume; 1104 } 1105 } 1106 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 1107 } 1108 if (useFVM) { 1109 ierr = VecRestoreArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1110 ierr = VecRestoreArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);CHKERRQ(ierr); 1111 ierr = VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1112 ierr = DMRestoreLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1113 ierr = VecDestroy(&faceGeometryFVM);CHKERRQ(ierr); 1114 ierr = VecDestroy(&cellGeometryFVM);CHKERRQ(ierr); 1115 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr); 1116 } 1117 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1118 if (mesh->printFEM) { 1119 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");CHKERRQ(ierr); 1120 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);CHKERRQ(ierr);} 1121 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 1122 } 1123 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1124 ierr = MPIU_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr); 1125 ierr = PetscFree3(lintegral,u,cgeomFEM);CHKERRQ(ierr); 1126 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1127 PetscFunctionReturn(0); 1128 } 1129 1130 /*@ 1131 DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1132 1133 Input Parameters: 1134 + dmf - The fine mesh 1135 . dmc - The coarse mesh 1136 - user - The user context 1137 1138 Output Parameter: 1139 . In - The interpolation matrix 1140 1141 Level: developer 1142 1143 .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1144 @*/ 1145 PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user) 1146 { 1147 DM_Plex *mesh = (DM_Plex *) dmc->data; 1148 const char *name = "Interpolator"; 1149 PetscDS prob; 1150 PetscFE *feRef; 1151 PetscFV *fvRef; 1152 PetscSection fsection, fglobalSection; 1153 PetscSection csection, cglobalSection; 1154 PetscScalar *elemMat; 1155 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c; 1156 PetscInt cTotDim, rTotDim = 0; 1157 PetscErrorCode ierr; 1158 1159 PetscFunctionBegin; 1160 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1161 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1162 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1163 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1164 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1165 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1166 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1167 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1168 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1169 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1170 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr); 1171 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1172 for (f = 0; f < Nf; ++f) { 1173 PetscObject obj; 1174 PetscClassId id; 1175 PetscInt rNb = 0, Nc = 0; 1176 1177 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1178 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1179 if (id == PETSCFE_CLASSID) { 1180 PetscFE fe = (PetscFE) obj; 1181 1182 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1183 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1184 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1185 } else if (id == PETSCFV_CLASSID) { 1186 PetscFV fv = (PetscFV) obj; 1187 PetscDualSpace Q; 1188 1189 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1190 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1191 ierr = PetscDualSpaceGetDimension(Q, &rNb);CHKERRQ(ierr); 1192 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1193 } 1194 rTotDim += rNb; 1195 } 1196 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1197 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 1198 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 1199 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1200 PetscDualSpace Qref; 1201 PetscQuadrature f; 1202 const PetscReal *qpoints, *qweights; 1203 PetscReal *points; 1204 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1205 1206 /* Compose points from all dual basis functionals */ 1207 if (feRef[fieldI]) { 1208 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1209 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 1210 } else { 1211 ierr = PetscFVGetDualSpace(fvRef[fieldI], &Qref);CHKERRQ(ierr); 1212 ierr = PetscFVGetNumComponents(fvRef[fieldI], &Nc);CHKERRQ(ierr); 1213 } 1214 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1215 for (i = 0; i < fpdim; ++i) { 1216 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1217 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1218 npoints += Np; 1219 } 1220 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1221 for (i = 0, k = 0; i < fpdim; ++i) { 1222 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1223 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1224 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1225 } 1226 1227 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1228 PetscObject obj; 1229 PetscClassId id; 1230 PetscReal *B; 1231 PetscInt NcJ = 0, cpdim = 0, j, qNc; 1232 1233 ierr = PetscDSGetDiscretization(prob, fieldJ, &obj);CHKERRQ(ierr); 1234 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1235 if (id == PETSCFE_CLASSID) { 1236 PetscFE fe = (PetscFE) obj; 1237 1238 /* Evaluate basis at points */ 1239 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 1240 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1241 /* For now, fields only interpolate themselves */ 1242 if (fieldI == fieldJ) { 1243 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); 1244 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1245 for (i = 0, k = 0; i < fpdim; ++i) { 1246 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1247 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);CHKERRQ(ierr); 1248 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); 1249 for (p = 0; p < Np; ++p, ++k) { 1250 for (j = 0; j < cpdim; ++j) { 1251 /* NOTE: This is not quite right, unless fpdim == number of fine grid functional quad points */ 1252 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p*qNc+c]; 1253 } 1254 } 1255 } 1256 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1257 } 1258 } else if (id == PETSCFV_CLASSID) { 1259 PetscFV fv = (PetscFV) obj; 1260 1261 /* Evaluate constant function at points */ 1262 ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr); 1263 cpdim = 1; 1264 /* For now, fields only interpolate themselves */ 1265 if (fieldI == fieldJ) { 1266 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); 1267 for (i = 0, k = 0; i < fpdim; ++i) { 1268 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1269 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);CHKERRQ(ierr); 1270 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); 1271 for (p = 0; p < Np; ++p, ++k) { 1272 for (j = 0; j < cpdim; ++j) { 1273 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p*qNc+c]; 1274 } 1275 } 1276 } 1277 } 1278 } 1279 offsetJ += cpdim*NcJ; 1280 } 1281 offsetI += fpdim*Nc; 1282 ierr = PetscFree(points);CHKERRQ(ierr); 1283 } 1284 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1285 /* Preallocate matrix */ 1286 { 1287 Mat preallocator; 1288 PetscScalar *vals; 1289 PetscInt *cellCIndices, *cellFIndices; 1290 PetscInt locRows, locCols, cell; 1291 1292 ierr = MatGetLocalSize(In, &locRows, &locCols);CHKERRQ(ierr); 1293 ierr = MatCreate(PetscObjectComm((PetscObject) In), &preallocator);CHKERRQ(ierr); 1294 ierr = MatSetType(preallocator, MATPREALLOCATOR);CHKERRQ(ierr); 1295 ierr = MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1296 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 1297 ierr = PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1298 for (cell = cStart; cell < cEnd; ++cell) { 1299 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1300 ierr = MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);CHKERRQ(ierr); 1301 } 1302 ierr = PetscFree3(vals,cellCIndices,cellFIndices);CHKERRQ(ierr); 1303 ierr = MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1304 ierr = MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1305 ierr = MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);CHKERRQ(ierr); 1306 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 1307 } 1308 /* Fill matrix */ 1309 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1310 for (c = cStart; c < cEnd; ++c) { 1311 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1312 } 1313 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1314 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1315 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1316 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1317 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1318 if (mesh->printFEM) { 1319 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1320 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1321 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1322 } 1323 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1324 PetscFunctionReturn(0); 1325 } 1326 1327 /*@ 1328 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 1329 1330 Input Parameters: 1331 + dmf - The fine mesh 1332 . dmc - The coarse mesh 1333 - user - The user context 1334 1335 Output Parameter: 1336 . In - The interpolation matrix 1337 1338 Level: developer 1339 1340 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1341 @*/ 1342 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 1343 { 1344 DM_Plex *mesh = (DM_Plex *) dmf->data; 1345 const char *name = "Interpolator"; 1346 PetscDS prob; 1347 PetscSection fsection, csection, globalFSection, globalCSection; 1348 PetscHashJK ht; 1349 PetscLayout rLayout; 1350 PetscInt *dnz, *onz; 1351 PetscInt locRows, rStart, rEnd; 1352 PetscReal *x, *v0, *J, *invJ, detJ; 1353 PetscReal *v0c, *Jc, *invJc, detJc; 1354 PetscScalar *elemMat; 1355 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1356 PetscErrorCode ierr; 1357 1358 PetscFunctionBegin; 1359 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1360 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1361 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1362 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1363 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1364 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1365 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1366 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1367 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1368 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1369 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1370 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1371 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1372 ierr = PetscMalloc1(totDim, &elemMat);CHKERRQ(ierr); 1373 1374 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1375 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1376 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1377 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1378 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1379 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1380 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1381 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1382 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1383 for (field = 0; field < Nf; ++field) { 1384 PetscObject obj; 1385 PetscClassId id; 1386 PetscDualSpace Q = NULL; 1387 PetscQuadrature f; 1388 const PetscReal *qpoints; 1389 PetscInt Nc, Np, fpdim, i, d; 1390 1391 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1392 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1393 if (id == PETSCFE_CLASSID) { 1394 PetscFE fe = (PetscFE) obj; 1395 1396 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1397 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1398 } else if (id == PETSCFV_CLASSID) { 1399 PetscFV fv = (PetscFV) obj; 1400 1401 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1402 Nc = 1; 1403 } 1404 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1405 /* For each fine grid cell */ 1406 for (cell = cStart; cell < cEnd; ++cell) { 1407 PetscInt *findices, *cindices; 1408 PetscInt numFIndices, numCIndices; 1409 1410 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1411 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1412 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim); 1413 for (i = 0; i < fpdim; ++i) { 1414 Vec pointVec; 1415 PetscScalar *pV; 1416 PetscSF coarseCellSF = NULL; 1417 const PetscSFNode *coarseCells; 1418 PetscInt numCoarseCells, q, c; 1419 1420 /* Get points from the dual basis functional quadrature */ 1421 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1422 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1423 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1424 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1425 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1426 for (q = 0; q < Np; ++q) { 1427 /* Transform point to real space */ 1428 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1429 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1430 } 1431 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1432 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1433 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1434 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 1435 /* Update preallocation info */ 1436 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1437 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1438 { 1439 PetscHashJKKey key; 1440 PetscHashJKIter missing, iter; 1441 1442 key.j = findices[i]; 1443 if (key.j >= 0) { 1444 /* Get indices for coarse elements */ 1445 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1446 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1447 for (c = 0; c < numCIndices; ++c) { 1448 key.k = cindices[c]; 1449 if (key.k < 0) continue; 1450 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1451 if (missing) { 1452 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1453 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1454 else ++onz[key.j-rStart]; 1455 } 1456 } 1457 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1458 } 1459 } 1460 } 1461 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1462 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1463 } 1464 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1465 } 1466 } 1467 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1468 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1469 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1470 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1471 for (field = 0; field < Nf; ++field) { 1472 PetscObject obj; 1473 PetscClassId id; 1474 PetscDualSpace Q = NULL; 1475 PetscQuadrature f; 1476 const PetscReal *qpoints, *qweights; 1477 PetscInt Nc, qNc, Np, fpdim, i, d; 1478 1479 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1480 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1481 if (id == PETSCFE_CLASSID) { 1482 PetscFE fe = (PetscFE) obj; 1483 1484 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1485 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1486 } else if (id == PETSCFV_CLASSID) { 1487 PetscFV fv = (PetscFV) obj; 1488 1489 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1490 Nc = 1; 1491 } 1492 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1493 /* For each fine grid cell */ 1494 for (cell = cStart; cell < cEnd; ++cell) { 1495 PetscInt *findices, *cindices; 1496 PetscInt numFIndices, numCIndices; 1497 1498 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1499 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1500 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim); 1501 for (i = 0; i < fpdim; ++i) { 1502 Vec pointVec; 1503 PetscScalar *pV; 1504 PetscSF coarseCellSF = NULL; 1505 const PetscSFNode *coarseCells; 1506 PetscInt numCoarseCells, cpdim, q, c, j; 1507 1508 /* Get points from the dual basis functional quadrature */ 1509 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1510 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);CHKERRQ(ierr); 1511 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); 1512 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1513 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1514 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1515 for (q = 0; q < Np; ++q) { 1516 /* Transform point to real space */ 1517 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1518 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1519 } 1520 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1521 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1522 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1523 /* Update preallocation info */ 1524 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1525 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1526 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1527 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1528 PetscReal pVReal[3]; 1529 1530 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1531 /* Transform points from real space to coarse reference space */ 1532 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 1533 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 1534 CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x); 1535 1536 if (id == PETSCFE_CLASSID) { 1537 PetscFE fe = (PetscFE) obj; 1538 PetscReal *B; 1539 1540 /* Evaluate coarse basis on contained point */ 1541 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1542 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 1543 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 1544 /* Get elemMat entries by multiplying by weight */ 1545 for (j = 0; j < cpdim; ++j) { 1546 for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*qweights[ccell*qNc + c]; 1547 } 1548 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1549 } else { 1550 cpdim = 1; 1551 for (j = 0; j < cpdim; ++j) { 1552 for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c]; 1553 } 1554 } 1555 /* Update interpolator */ 1556 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 1557 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 1558 ierr = MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1559 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1560 } 1561 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1562 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1563 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1564 } 1565 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1566 } 1567 } 1568 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 1569 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 1570 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1571 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1572 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1573 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1574 PetscFunctionReturn(0); 1575 } 1576 1577 /*@ 1578 DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns 1579 1580 Input Parameters: 1581 + dmc - The coarse mesh 1582 - dmf - The fine mesh 1583 - user - The user context 1584 1585 Output Parameter: 1586 . sc - The mapping 1587 1588 Level: developer 1589 1590 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1591 @*/ 1592 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 1593 { 1594 PetscDS prob; 1595 PetscFE *feRef; 1596 PetscFV *fvRef; 1597 Vec fv, cv; 1598 IS fis, cis; 1599 PetscSection fsection, fglobalSection, csection, cglobalSection; 1600 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 1601 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m; 1602 PetscErrorCode ierr; 1603 1604 PetscFunctionBegin; 1605 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1606 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1607 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1608 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1609 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1610 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1611 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1612 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1613 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1614 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1615 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1616 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1617 for (f = 0; f < Nf; ++f) { 1618 PetscObject obj; 1619 PetscClassId id; 1620 PetscInt fNb = 0, Nc = 0; 1621 1622 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1623 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1624 if (id == PETSCFE_CLASSID) { 1625 PetscFE fe = (PetscFE) obj; 1626 1627 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1628 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 1629 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1630 } else if (id == PETSCFV_CLASSID) { 1631 PetscFV fv = (PetscFV) obj; 1632 PetscDualSpace Q; 1633 1634 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1635 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1636 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 1637 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1638 } 1639 fTotDim += fNb*Nc; 1640 } 1641 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1642 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 1643 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 1644 PetscFE feC; 1645 PetscFV fvC; 1646 PetscDualSpace QF, QC; 1647 PetscInt NcF, NcC, fpdim, cpdim; 1648 1649 if (feRef[field]) { 1650 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 1651 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 1652 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 1653 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 1654 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1655 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 1656 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1657 } else { 1658 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 1659 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 1660 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 1661 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 1662 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1663 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 1664 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1665 } 1666 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); 1667 for (c = 0; c < cpdim; ++c) { 1668 PetscQuadrature cfunc; 1669 const PetscReal *cqpoints; 1670 PetscInt NpC; 1671 PetscBool found = PETSC_FALSE; 1672 1673 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 1674 ierr = PetscQuadratureGetData(cfunc, NULL, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 1675 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 1676 for (f = 0; f < fpdim; ++f) { 1677 PetscQuadrature ffunc; 1678 const PetscReal *fqpoints; 1679 PetscReal sum = 0.0; 1680 PetscInt NpF, comp; 1681 1682 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 1683 ierr = PetscQuadratureGetData(ffunc, NULL, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 1684 if (NpC != NpF) continue; 1685 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 1686 if (sum > 1.0e-9) continue; 1687 for (comp = 0; comp < NcC; ++comp) { 1688 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 1689 } 1690 found = PETSC_TRUE; 1691 break; 1692 } 1693 if (!found) { 1694 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 1695 if (fvRef[field]) { 1696 PetscInt comp; 1697 for (comp = 0; comp < NcC; ++comp) { 1698 cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp; 1699 } 1700 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 1701 } 1702 } 1703 offsetC += cpdim*NcC; 1704 offsetF += fpdim*NcF; 1705 } 1706 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 1707 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1708 1709 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 1710 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 1711 ierr = VecGetOwnershipRange(cv, &startC, &endC);CHKERRQ(ierr); 1712 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 1713 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 1714 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 1715 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 1716 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 1717 for (c = cStart; c < cEnd; ++c) { 1718 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 1719 for (d = 0; d < cTotDim; ++d) { 1720 if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue; 1721 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]]); 1722 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 1723 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 1724 } 1725 } 1726 ierr = PetscFree(cmap);CHKERRQ(ierr); 1727 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 1728 1729 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 1730 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 1731 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 1732 ierr = ISDestroy(&cis);CHKERRQ(ierr); 1733 ierr = ISDestroy(&fis);CHKERRQ(ierr); 1734 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 1735 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 1736 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1737 PetscFunctionReturn(0); 1738 } 1739