1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petscsf.h> 3 4 #include <petsc/private/petscfeimpl.h> 5 #include <petsc/private/petscfvimpl.h> 6 7 /*@ 8 DMPlexGetScale - Get the scale for the specified fundamental unit 9 10 Not collective 11 12 Input Arguments: 13 + dm - the DM 14 - unit - The SI unit 15 16 Output Argument: 17 . scale - The value used to scale all quantities with this unit 18 19 Level: advanced 20 21 .seealso: DMPlexSetScale(), PetscUnit 22 @*/ 23 PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale) 24 { 25 DM_Plex *mesh = (DM_Plex*) dm->data; 26 27 PetscFunctionBegin; 28 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 29 PetscValidPointer(scale, 3); 30 *scale = mesh->scale[unit]; 31 PetscFunctionReturn(0); 32 } 33 34 /*@ 35 DMPlexSetScale - Set the scale for the specified fundamental unit 36 37 Not collective 38 39 Input Arguments: 40 + dm - the DM 41 . unit - The SI unit 42 - scale - The value used to scale all quantities with this unit 43 44 Level: advanced 45 46 .seealso: DMPlexGetScale(), PetscUnit 47 @*/ 48 PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale) 49 { 50 DM_Plex *mesh = (DM_Plex*) dm->data; 51 52 PetscFunctionBegin; 53 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 54 mesh->scale[unit] = scale; 55 PetscFunctionReturn(0); 56 } 57 58 PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k) 59 { 60 switch (i) { 61 case 0: 62 switch (j) { 63 case 0: return 0; 64 case 1: 65 switch (k) { 66 case 0: return 0; 67 case 1: return 0; 68 case 2: return 1; 69 } 70 case 2: 71 switch (k) { 72 case 0: return 0; 73 case 1: return -1; 74 case 2: return 0; 75 } 76 } 77 case 1: 78 switch (j) { 79 case 0: 80 switch (k) { 81 case 0: return 0; 82 case 1: return 0; 83 case 2: return -1; 84 } 85 case 1: return 0; 86 case 2: 87 switch (k) { 88 case 0: return 1; 89 case 1: return 0; 90 case 2: return 0; 91 } 92 } 93 case 2: 94 switch (j) { 95 case 0: 96 switch (k) { 97 case 0: return 0; 98 case 1: return 1; 99 case 2: return 0; 100 } 101 case 1: 102 switch (k) { 103 case 0: return -1; 104 case 1: return 0; 105 case 2: return 0; 106 } 107 case 2: return 0; 108 } 109 } 110 return 0; 111 } 112 113 static PetscErrorCode DMPlexProjectRigidBody_Private(PetscInt dim, PetscReal t, const PetscReal X[], PetscInt Nf, PetscScalar *mode, void *ctx) 114 { 115 PetscInt *ctxInt = (PetscInt *) ctx; 116 PetscInt dim2 = ctxInt[0]; 117 PetscInt d = ctxInt[1]; 118 PetscInt i, j, k = dim > 2 ? d - dim : d; 119 120 PetscFunctionBegin; 121 if (dim != dim2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input dimension %d does not match context dimension %d", dim, dim2); 122 for (i = 0; i < dim; i++) mode[i] = 0.; 123 if (d < dim) { 124 mode[d] = 1.; 125 } else { 126 for (i = 0; i < dim; i++) { 127 for (j = 0; j < dim; j++) { 128 mode[j] += epsilon(i, j, k)*X[i]; 129 } 130 } 131 } 132 PetscFunctionReturn(0); 133 } 134 135 /*@C 136 DMPlexCreateRigidBody - for the default global section, create rigid body modes from coordinates 137 138 Collective on DM 139 140 Input Arguments: 141 . dm - the DM 142 143 Output Argument: 144 . sp - the null space 145 146 Note: This is necessary to take account of Dirichlet conditions on the displacements 147 148 Level: advanced 149 150 .seealso: MatNullSpaceCreate() 151 @*/ 152 PetscErrorCode DMPlexCreateRigidBody(DM dm, MatNullSpace *sp) 153 { 154 MPI_Comm comm; 155 Vec mode[6]; 156 PetscSection section, globalSection; 157 PetscInt dim, dimEmbed, n, m, d, i, j; 158 PetscErrorCode ierr; 159 160 PetscFunctionBegin; 161 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 162 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 163 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 164 if (dim == 1) { 165 ierr = MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);CHKERRQ(ierr); 166 PetscFunctionReturn(0); 167 } 168 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 169 ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr); 170 ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr); 171 m = (dim*(dim+1))/2; 172 ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr); 173 ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr); 174 ierr = VecSetUp(mode[0]);CHKERRQ(ierr); 175 for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);} 176 for (d = 0; d < m; d++) { 177 PetscInt ctx[2]; 178 PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private; 179 void *voidctx = (void *) (&ctx[0]); 180 181 ctx[0] = dimEmbed; 182 ctx[1] = d; 183 ierr = DMProjectFunction(dm, 0.0, &func, &voidctx, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 184 } 185 for (i = 0; i < dim; ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);} 186 /* Orthonormalize system */ 187 for (i = dim; i < m; ++i) { 188 PetscScalar dots[6]; 189 190 ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr); 191 for (j = 0; j < i; ++j) dots[j] *= -1.0; 192 ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr); 193 ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr); 194 } 195 ierr = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);CHKERRQ(ierr); 196 for (i = 0; i< m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);} 197 PetscFunctionReturn(0); 198 } 199 200 /*@ 201 DMPlexSetMaxProjectionHeight - In DMPlexProjectXXXLocal() functions, the projected values of a basis function's dofs 202 are computed by associating the basis function with one of the mesh points in its transitively-closed support, and 203 evaluating the dual space basis of that point. A basis function is associated with the point in its 204 transitively-closed support whose mesh height is highest (w.r.t. DAG height), but not greater than the maximum 205 projection height, which is set with this function. By default, the maximum projection height is zero, which means 206 that only mesh cells are used to project basis functions. A height of one, for example, evaluates a cell-interior 207 basis functions using its cells dual space basis, but all other basis functions with the dual space basis of a face. 208 209 Input Parameters: 210 + dm - the DMPlex object 211 - height - the maximum projection height >= 0 212 213 Level: advanced 214 215 .seealso: DMPlexGetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal() 216 @*/ 217 PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height) 218 { 219 DM_Plex *plex = (DM_Plex *) dm->data; 220 221 PetscFunctionBegin; 222 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 223 plex->maxProjectionHeight = height; 224 PetscFunctionReturn(0); 225 } 226 227 /*@ 228 DMPlexGetMaxProjectionHeight - Get the maximum height (w.r.t. DAG) of mesh points used to evaluate dual bases in 229 DMPlexProjectXXXLocal() functions. 230 231 Input Parameters: 232 . dm - the DMPlex object 233 234 Output Parameters: 235 . height - the maximum projection height 236 237 Level: intermediate 238 239 .seealso: DMPlexSetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal() 240 @*/ 241 PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height) 242 { 243 DM_Plex *plex = (DM_Plex *) dm->data; 244 245 PetscFunctionBegin; 246 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 247 *height = plex->maxProjectionHeight; 248 PetscFunctionReturn(0); 249 } 250 251 /*@C 252 DMPlexInsertBoundaryValuesEssential - Insert boundary values into a local vector 253 254 Input Parameters: 255 + dm - The DM, with a PetscDS that matches the problem being constrained 256 . time - The time 257 . field - The field to constrain 258 . Nc - The number of constrained field components, or 0 for all components 259 . comps - An array of constrained component numbers, or NULL for all components 260 . label - The DMLabel defining constrained points 261 . numids - The number of DMLabel ids for constrained points 262 . ids - An array of ids for constrained points 263 . func - A pointwise function giving boundary values 264 - ctx - An optional user context for bcFunc 265 266 Output Parameter: 267 . locX - A local vector to receives the boundary values 268 269 Level: developer 270 271 .seealso: DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary() 272 @*/ 273 PetscErrorCode DMPlexInsertBoundaryValuesEssential(DM dm, PetscReal time, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void *ctx, Vec locX) 274 { 275 PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx); 276 void **ctxs; 277 PetscInt numFields; 278 PetscErrorCode ierr; 279 280 PetscFunctionBegin; 281 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 282 ierr = PetscCalloc2(numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 283 funcs[field] = func; 284 ctxs[field] = ctx; 285 ierr = DMProjectFunctionLabelLocal(dm, time, label, numids, ids, Nc, comps, funcs, ctxs, INSERT_BC_VALUES, locX);CHKERRQ(ierr); 286 ierr = PetscFree2(funcs,ctxs);CHKERRQ(ierr); 287 PetscFunctionReturn(0); 288 } 289 290 /*@C 291 DMPlexInsertBoundaryValuesEssentialField - Insert boundary values into a local vector 292 293 Input Parameters: 294 + dm - The DM, with a PetscDS that matches the problem being constrained 295 . time - The time 296 . locU - A local vector with the input solution values 297 . field - The field to constrain 298 . Nc - The number of constrained field components, or 0 for all components 299 . comps - An array of constrained component numbers, or NULL for all components 300 . label - The DMLabel defining constrained points 301 . numids - The number of DMLabel ids for constrained points 302 . ids - An array of ids for constrained points 303 . func - A pointwise function giving boundary values 304 - ctx - An optional user context for bcFunc 305 306 Output Parameter: 307 . locX - A local vector to receives the boundary values 308 309 Level: developer 310 311 .seealso: DMPlexInsertBoundaryValuesEssential(), DMAddBoundary() 312 @*/ 313 PetscErrorCode DMPlexInsertBoundaryValuesEssentialField(DM dm, PetscReal time, Vec locU, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[], 314 void (*func)(PetscInt, PetscInt, PetscInt, 315 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 316 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 317 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], 318 PetscScalar[]), 319 void *ctx, Vec locX) 320 { 321 void (**funcs)(PetscInt, PetscInt, PetscInt, 322 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 323 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 324 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]); 325 void **ctxs; 326 PetscInt numFields; 327 PetscErrorCode ierr; 328 329 PetscFunctionBegin; 330 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 331 ierr = PetscCalloc2(numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 332 funcs[field] = func; 333 ctxs[field] = ctx; 334 ierr = DMProjectFieldLabelLocal(dm, time, label, numids, ids, Nc, comps, locU, funcs, INSERT_BC_VALUES, locX);CHKERRQ(ierr); 335 ierr = PetscFree2(funcs,ctxs);CHKERRQ(ierr); 336 PetscFunctionReturn(0); 337 } 338 339 /*@C 340 DMPlexInsertBoundaryValuesRiemann - Insert boundary values into a local vector 341 342 Input Parameters: 343 + dm - The DM, with a PetscDS that matches the problem being constrained 344 . time - The time 345 . faceGeometry - A vector with the FVM face geometry information 346 . cellGeometry - A vector with the FVM cell geometry information 347 . Grad - A vector with the FVM cell gradient information 348 . field - The field to constrain 349 . Nc - The number of constrained field components, or 0 for all components 350 . comps - An array of constrained component numbers, or NULL for all components 351 . label - The DMLabel defining constrained points 352 . numids - The number of DMLabel ids for constrained points 353 . ids - An array of ids for constrained points 354 . func - A pointwise function giving boundary values 355 - ctx - An optional user context for bcFunc 356 357 Output Parameter: 358 . locX - A local vector to receives the boundary values 359 360 Note: This implementation currently ignores the numcomps/comps argument from DMAddBoundary() 361 362 Level: developer 363 364 .seealso: DMPlexInsertBoundaryValuesEssential(), DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary() 365 @*/ 366 PetscErrorCode DMPlexInsertBoundaryValuesRiemann(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[], 367 PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX) 368 { 369 PetscDS prob; 370 PetscSF sf; 371 DM dmFace, dmCell, dmGrad; 372 const PetscScalar *facegeom, *cellgeom = NULL, *grad; 373 const PetscInt *leaves; 374 PetscScalar *x, *fx; 375 PetscInt dim, nleaves, loc, fStart, fEnd, pdim, i; 376 PetscErrorCode ierr, ierru = 0; 377 378 PetscFunctionBegin; 379 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 380 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 381 nleaves = PetscMax(0, nleaves); 382 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 383 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 384 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 385 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 386 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 387 if (cellGeometry) { 388 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 389 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 390 } 391 if (Grad) { 392 PetscFV fv; 393 394 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);CHKERRQ(ierr); 395 ierr = VecGetDM(Grad, &dmGrad);CHKERRQ(ierr); 396 ierr = VecGetArrayRead(Grad, &grad);CHKERRQ(ierr); 397 ierr = PetscFVGetNumComponents(fv, &pdim);CHKERRQ(ierr); 398 ierr = DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 399 } 400 ierr = VecGetArray(locX, &x);CHKERRQ(ierr); 401 for (i = 0; i < numids; ++i) { 402 IS faceIS; 403 const PetscInt *faces; 404 PetscInt numFaces, f; 405 406 ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr); 407 if (!faceIS) continue; /* No points with that id on this process */ 408 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 409 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 410 for (f = 0; f < numFaces; ++f) { 411 const PetscInt face = faces[f], *cells; 412 PetscFVFaceGeom *fg; 413 414 if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */ 415 ierr = PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);CHKERRQ(ierr); 416 if (loc >= 0) continue; 417 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 418 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 419 if (Grad) { 420 PetscFVCellGeom *cg; 421 PetscScalar *cx, *cgrad; 422 PetscScalar *xG; 423 PetscReal dx[3]; 424 PetscInt d; 425 426 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);CHKERRQ(ierr); 427 ierr = DMPlexPointLocalRead(dm, cells[0], x, &cx);CHKERRQ(ierr); 428 ierr = DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);CHKERRQ(ierr); 429 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 430 DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx); 431 for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx); 432 ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx); 433 if (ierru) { 434 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 435 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 436 goto cleanup; 437 } 438 } else { 439 PetscScalar *xI; 440 PetscScalar *xG; 441 442 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr); 443 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 444 ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx); 445 if (ierru) { 446 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 447 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 448 goto cleanup; 449 } 450 } 451 } 452 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 453 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 454 } 455 cleanup: 456 ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr); 457 if (Grad) { 458 ierr = DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 459 ierr = VecRestoreArrayRead(Grad, &grad);CHKERRQ(ierr); 460 } 461 if (cellGeometry) {ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr);} 462 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 463 CHKERRQ(ierru); 464 PetscFunctionReturn(0); 465 } 466 467 PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM) 468 { 469 PetscInt numBd, b; 470 PetscErrorCode ierr; 471 472 PetscFunctionBegin; 473 ierr = PetscDSGetNumBoundary(dm->prob, &numBd);CHKERRQ(ierr); 474 for (b = 0; b < numBd; ++b) { 475 DMBoundaryConditionType type; 476 const char *labelname; 477 DMLabel label; 478 PetscInt field, Nc; 479 const PetscInt *comps; 480 PetscObject obj; 481 PetscClassId id; 482 void (*func)(void); 483 PetscInt numids; 484 const PetscInt *ids; 485 void *ctx; 486 487 ierr = DMGetBoundary(dm, b, &type, NULL, &labelname, &field, &Nc, &comps, &func, &numids, &ids, &ctx);CHKERRQ(ierr); 488 if (insertEssential != (type & DM_BC_ESSENTIAL)) continue; 489 ierr = DMGetLabel(dm, labelname, &label);CHKERRQ(ierr); 490 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 491 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 492 if (id == PETSCFE_CLASSID) { 493 switch (type) { 494 /* for FEM, there is no insertion to be done for non-essential boundary conditions */ 495 case DM_BC_ESSENTIAL: 496 ierr = DMPlexLabelAddCells(dm,label);CHKERRQ(ierr); 497 ierr = DMPlexInsertBoundaryValuesEssential(dm, time, field, Nc, comps, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);CHKERRQ(ierr); 498 ierr = DMPlexLabelClearCells(dm,label);CHKERRQ(ierr); 499 break; 500 case DM_BC_ESSENTIAL_FIELD: 501 ierr = DMPlexLabelAddCells(dm,label);CHKERRQ(ierr); 502 ierr = DMPlexInsertBoundaryValuesEssentialField(dm, time, locX, field, Nc, comps, label, numids, ids, 503 (void (*)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 504 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 505 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) func, ctx, locX);CHKERRQ(ierr); 506 ierr = DMPlexLabelClearCells(dm,label);CHKERRQ(ierr); 507 break; 508 default: break; 509 } 510 } else if (id == PETSCFV_CLASSID) { 511 if (!faceGeomFVM) continue; 512 ierr = DMPlexInsertBoundaryValuesRiemann(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, field, Nc, comps, label, numids, ids, 513 (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);CHKERRQ(ierr); 514 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 515 } 516 PetscFunctionReturn(0); 517 } 518 519 /*@ 520 DMPlexInsertBoundaryValues - Puts coefficients which represent boundary values into the local solution vector 521 522 Input Parameters: 523 + dm - The DM 524 . insertEssential - Should I insert essential (e.g. Dirichlet) or inessential (e.g. Neumann) boundary conditions 525 . time - The time 526 . faceGeomFVM - Face geometry data for FV discretizations 527 . cellGeomFVM - Cell geometry data for FV discretizations 528 - gradFVM - Gradient reconstruction data for FV discretizations 529 530 Output Parameters: 531 . locX - Solution updated with boundary values 532 533 Level: developer 534 535 .seealso: DMProjectFunctionLabelLocal() 536 @*/ 537 PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM) 538 { 539 PetscErrorCode ierr; 540 541 PetscFunctionBegin; 542 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 543 PetscValidHeaderSpecific(locX, VEC_CLASSID, 2); 544 if (faceGeomFVM) {PetscValidHeaderSpecific(faceGeomFVM, VEC_CLASSID, 4);} 545 if (cellGeomFVM) {PetscValidHeaderSpecific(cellGeomFVM, VEC_CLASSID, 5);} 546 if (gradFVM) {PetscValidHeaderSpecific(gradFVM, VEC_CLASSID, 6);} 547 ierr = PetscTryMethod(dm,"DMPlexInsertBoundaryValues_C",(DM,PetscBool,Vec,PetscReal,Vec,Vec,Vec),(dm,insertEssential,locX,time,faceGeomFVM,cellGeomFVM,gradFVM));CHKERRQ(ierr); 548 PetscFunctionReturn(0); 549 } 550 551 PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 552 { 553 const PetscInt debug = 0; 554 PetscSection section; 555 PetscQuadrature quad; 556 Vec localX; 557 PetscScalar *funcVal, *interpolant; 558 PetscReal *coords, *detJ, *J; 559 PetscReal localDiff = 0.0; 560 const PetscReal *quadWeights; 561 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset; 562 PetscErrorCode ierr; 563 564 PetscFunctionBegin; 565 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 566 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 567 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 568 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 569 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 570 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 571 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 572 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 573 for (field = 0; field < numFields; ++field) { 574 PetscObject obj; 575 PetscClassId id; 576 PetscInt Nc; 577 578 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 579 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 580 if (id == PETSCFE_CLASSID) { 581 PetscFE fe = (PetscFE) obj; 582 583 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 584 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 585 } else if (id == PETSCFV_CLASSID) { 586 PetscFV fv = (PetscFV) obj; 587 588 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 589 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 590 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 591 numComponents += Nc; 592 } 593 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);CHKERRQ(ierr); 594 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 595 ierr = PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);CHKERRQ(ierr); 596 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 597 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 598 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 599 for (c = cStart; c < cEnd; ++c) { 600 PetscScalar *x = NULL; 601 PetscReal elemDiff = 0.0; 602 PetscInt qc = 0; 603 604 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);CHKERRQ(ierr); 605 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 606 607 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 608 PetscObject obj; 609 PetscClassId id; 610 void * const ctx = ctxs ? ctxs[field] : NULL; 611 PetscInt Nb, Nc, q, fc; 612 613 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 614 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 615 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 616 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 617 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 618 if (debug) { 619 char title[1024]; 620 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 621 ierr = DMPrintCellVector(c, title, Nb, &x[fieldOffset]);CHKERRQ(ierr); 622 } 623 for (q = 0; q < Nq; ++q) { 624 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); 625 ierr = (*funcs[field])(coordDim, time, &coords[coordDim * q], Nc, funcVal, ctx); 626 if (ierr) { 627 PetscErrorCode ierr2; 628 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 629 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 630 ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2); 631 CHKERRQ(ierr); 632 } 633 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 634 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 635 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 636 for (fc = 0; fc < Nc; ++fc) { 637 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 638 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);} 639 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]; 640 } 641 } 642 fieldOffset += Nb; 643 qc += Nc; 644 } 645 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 646 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 647 localDiff += elemDiff; 648 } 649 ierr = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr); 650 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 651 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 652 *diff = PetscSqrtReal(*diff); 653 PetscFunctionReturn(0); 654 } 655 656 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) 657 { 658 const PetscInt debug = 0; 659 PetscSection section; 660 PetscQuadrature quad; 661 Vec localX; 662 PetscScalar *funcVal, *interpolant; 663 const PetscReal *quadPoints, *quadWeights; 664 PetscReal *coords, *realSpaceDer, *J, *invJ, *detJ; 665 PetscReal localDiff = 0.0; 666 PetscInt dim, coordDim, qNc = 0, Nq = 0, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset; 667 PetscErrorCode ierr; 668 669 PetscFunctionBegin; 670 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 671 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 672 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 673 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 674 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 675 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 676 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 677 for (field = 0; field < numFields; ++field) { 678 PetscFE fe; 679 PetscInt Nc; 680 681 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 682 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 683 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 684 numComponents += Nc; 685 } 686 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 687 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 688 /* ierr = DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 689 ierr = PetscMalloc7(numComponents,&funcVal,coordDim*Nq,&coords,coordDim*Nq,&realSpaceDer,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ,numComponents,&interpolant,Nq,&detJ);CHKERRQ(ierr); 690 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 691 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 692 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 693 for (c = cStart; c < cEnd; ++c) { 694 PetscScalar *x = NULL; 695 PetscReal elemDiff = 0.0; 696 PetscInt qc = 0; 697 698 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, invJ, detJ);CHKERRQ(ierr); 699 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 700 701 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 702 PetscFE fe; 703 void * const ctx = ctxs ? ctxs[field] : NULL; 704 PetscReal *basisDer; 705 PetscInt Nb, Nc, q, fc; 706 707 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 708 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 709 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 710 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 711 if (debug) { 712 char title[1024]; 713 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 714 ierr = DMPrintCellVector(c, title, Nb, &x[fieldOffset]);CHKERRQ(ierr); 715 } 716 for (q = 0; q < Nq; ++q) { 717 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); 718 ierr = (*funcs[field])(coordDim, time, &coords[q*coordDim], n, numFields, funcVal, ctx); 719 if (ierr) { 720 PetscErrorCode ierr2; 721 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 722 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 723 ierr2 = PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);CHKERRQ(ierr2); 724 CHKERRQ(ierr); 725 } 726 ierr = PetscFEInterpolateGradient_Static(fe, &x[fieldOffset], coordDim, invJ, n, q, interpolant);CHKERRQ(ierr); 727 for (fc = 0; fc < Nc; ++fc) { 728 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 729 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);} 730 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]; 731 } 732 } 733 fieldOffset += Nb; 734 qc += Nc; 735 } 736 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 737 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 738 localDiff += elemDiff; 739 } 740 ierr = PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);CHKERRQ(ierr); 741 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 742 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 743 *diff = PetscSqrtReal(*diff); 744 PetscFunctionReturn(0); 745 } 746 747 PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 748 { 749 const PetscInt debug = 0; 750 PetscSection section; 751 PetscQuadrature quad; 752 Vec localX; 753 PetscScalar *funcVal, *interpolant; 754 PetscReal *coords, *detJ, *J; 755 PetscReal *localDiff; 756 const PetscReal *quadPoints, *quadWeights; 757 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset; 758 PetscErrorCode ierr; 759 760 PetscFunctionBegin; 761 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 762 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 763 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 764 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 765 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 766 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 767 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 768 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 769 for (field = 0; field < numFields; ++field) { 770 PetscObject obj; 771 PetscClassId id; 772 PetscInt Nc; 773 774 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 775 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 776 if (id == PETSCFE_CLASSID) { 777 PetscFE fe = (PetscFE) obj; 778 779 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 780 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 781 } else if (id == PETSCFV_CLASSID) { 782 PetscFV fv = (PetscFV) obj; 783 784 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 785 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 786 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 787 numComponents += Nc; 788 } 789 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 790 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 791 ierr = PetscCalloc6(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);CHKERRQ(ierr); 792 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 793 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 794 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 795 for (c = cStart; c < cEnd; ++c) { 796 PetscScalar *x = NULL; 797 PetscInt qc = 0; 798 799 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);CHKERRQ(ierr); 800 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 801 802 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 803 PetscObject obj; 804 PetscClassId id; 805 void * const ctx = ctxs ? ctxs[field] : NULL; 806 PetscInt Nb, Nc, q, fc; 807 808 PetscReal elemDiff = 0.0; 809 810 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 811 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 812 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 813 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 814 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 815 if (debug) { 816 char title[1024]; 817 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 818 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 819 } 820 for (q = 0; q < Nq; ++q) { 821 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); 822 ierr = (*funcs[field])(coordDim, time, &coords[coordDim*q], numFields, funcVal, ctx); 823 if (ierr) { 824 PetscErrorCode ierr2; 825 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 826 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 827 ierr2 = PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2); 828 CHKERRQ(ierr); 829 } 830 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 831 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 832 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 833 for (fc = 0; fc < Nc; ++fc) { 834 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 835 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);} 836 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]; 837 } 838 } 839 fieldOffset += Nb; 840 qc += Nc; 841 localDiff[field] += elemDiff; 842 } 843 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 844 } 845 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 846 ierr = MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 847 for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]); 848 ierr = PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr); 849 PetscFunctionReturn(0); 850 } 851 852 /*@C 853 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. 854 855 Input Parameters: 856 + dm - The DM 857 . time - The time 858 . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation 859 . ctxs - Optional array of contexts to pass to each function, or NULL. 860 - X - The coefficient vector u_h 861 862 Output Parameter: 863 . D - A Vec which holds the difference ||u - u_h||_2 for each cell 864 865 Level: developer 866 867 .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff() 868 @*/ 869 PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D) 870 { 871 PetscSection section; 872 PetscQuadrature quad; 873 Vec localX; 874 PetscScalar *funcVal, *interpolant; 875 PetscReal *coords, *detJ, *J; 876 const PetscReal *quadPoints, *quadWeights; 877 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset; 878 PetscErrorCode ierr; 879 880 PetscFunctionBegin; 881 ierr = VecSet(D, 0.0);CHKERRQ(ierr); 882 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 883 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 884 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 885 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 886 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 887 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 888 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 889 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 890 for (field = 0; field < numFields; ++field) { 891 PetscObject obj; 892 PetscClassId id; 893 PetscInt Nc; 894 895 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 896 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 897 if (id == PETSCFE_CLASSID) { 898 PetscFE fe = (PetscFE) obj; 899 900 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 901 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 902 } else if (id == PETSCFV_CLASSID) { 903 PetscFV fv = (PetscFV) obj; 904 905 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 906 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 907 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 908 numComponents += Nc; 909 } 910 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 911 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 912 ierr = PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);CHKERRQ(ierr); 913 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 914 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 915 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 916 for (c = cStart; c < cEnd; ++c) { 917 PetscScalar *x = NULL; 918 PetscScalar elemDiff = 0.0; 919 PetscInt qc = 0; 920 921 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);CHKERRQ(ierr); 922 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 923 924 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 925 PetscObject obj; 926 PetscClassId id; 927 void * const ctx = ctxs ? ctxs[field] : NULL; 928 PetscInt Nb, Nc, q, fc; 929 930 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 931 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 932 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 933 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 934 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 935 if (funcs[field]) { 936 for (q = 0; q < Nq; ++q) { 937 if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", (double)detJ[q], c, q); 938 ierr = (*funcs[field])(dim, time, &coords[q*coordDim], Nc, funcVal, ctx); 939 if (ierr) { 940 PetscErrorCode ierr2; 941 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 942 ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2); 943 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 944 CHKERRQ(ierr); 945 } 946 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 947 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 948 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 949 for (fc = 0; fc < Nc; ++fc) { 950 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 951 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]; 952 } 953 } 954 } 955 fieldOffset += Nb; 956 qc += Nc; 957 } 958 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 959 ierr = VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);CHKERRQ(ierr); 960 } 961 ierr = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr); 962 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 963 ierr = VecSqrtAbs(D);CHKERRQ(ierr); 964 PetscFunctionReturn(0); 965 } 966 967 /*@C 968 DMPlexComputeGradientClementInterpolant - This function computes the L2 projection of the cellwise gradient of a function u onto P1, and stores it in a Vec. 969 970 Input Parameters: 971 + dm - The DM 972 - LocX - The coefficient vector u_h 973 974 Output Parameter: 975 . locC - A Vec which holds the Clement interpolant of the gradient 976 977 Notes: Add citation to (Clement, 1975) and definition of the interpolant 978 \nabla u_h(v_i) = \sum_{T_i \in support(v_i)} |T_i| \nabla u_h(T_i) / \sum_{T_i \in support(v_i)} |T_i| where |T_i| is the cell volume 979 980 Level: developer 981 982 .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff() 983 @*/ 984 PetscErrorCode DMPlexComputeGradientClementInterpolant(DM dm, Vec locX, Vec locC) 985 { 986 DM_Plex *mesh = (DM_Plex *) dm->data; 987 PetscInt debug = mesh->printFEM; 988 DM dmC; 989 PetscSection section; 990 PetscQuadrature quad; 991 PetscScalar *interpolant, *gradsum; 992 PetscReal *coords, *detJ, *J, *invJ; 993 const PetscReal *quadPoints, *quadWeights; 994 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, vStart, vEnd, v, field, fieldOffset; 995 PetscErrorCode ierr; 996 997 PetscFunctionBegin; 998 ierr = VecGetDM(locC, &dmC);CHKERRQ(ierr); 999 ierr = VecSet(locC, 0.0);CHKERRQ(ierr); 1000 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1001 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 1002 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1003 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1004 for (field = 0; field < numFields; ++field) { 1005 PetscObject obj; 1006 PetscClassId id; 1007 PetscInt Nc; 1008 1009 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1010 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1011 if (id == PETSCFE_CLASSID) { 1012 PetscFE fe = (PetscFE) obj; 1013 1014 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1015 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1016 } else if (id == PETSCFV_CLASSID) { 1017 PetscFV fv = (PetscFV) obj; 1018 1019 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1020 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1021 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1022 numComponents += Nc; 1023 } 1024 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 1025 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 1026 ierr = PetscMalloc6(coordDim*numComponents*2,&gradsum,coordDim*numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ);CHKERRQ(ierr); 1027 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1028 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1029 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1030 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1031 for (v = vStart; v < vEnd; ++v) { 1032 PetscScalar volsum = 0.0; 1033 PetscInt *star = NULL; 1034 PetscInt starSize, st, d, fc; 1035 1036 ierr = PetscMemzero(gradsum, coordDim*numComponents * sizeof(PetscScalar));CHKERRQ(ierr); 1037 ierr = DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1038 for (st = 0; st < starSize*2; st += 2) { 1039 const PetscInt cell = star[st]; 1040 PetscScalar *grad = &gradsum[coordDim*numComponents]; 1041 PetscScalar *x = NULL; 1042 PetscReal vol = 0.0; 1043 1044 if ((cell < cStart) || (cell >= cEnd)) continue; 1045 ierr = DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, J, invJ, detJ);CHKERRQ(ierr); 1046 ierr = DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr); 1047 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1048 PetscObject obj; 1049 PetscClassId id; 1050 PetscInt Nb, Nc, q, qc = 0; 1051 1052 ierr = PetscMemzero(grad, coordDim*numComponents * sizeof(PetscScalar));CHKERRQ(ierr); 1053 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1054 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1055 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1056 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1057 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1058 for (q = 0; q < Nq; ++q) { 1059 if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", (double)detJ[q], cell, q); 1060 if (ierr) { 1061 PetscErrorCode ierr2; 1062 ierr2 = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr2); 1063 ierr2 = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr2); 1064 ierr2 = PetscFree4(interpolant,coords,detJ,J);CHKERRQ(ierr2); 1065 CHKERRQ(ierr); 1066 } 1067 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolateGradient_Static((PetscFE) obj, &x[fieldOffset], coordDim, invJ, NULL, q, interpolant);CHKERRQ(ierr);} 1068 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1069 for (fc = 0; fc < Nc; ++fc) { 1070 const PetscReal wt = quadWeights[q*qNc+qc+fc]; 1071 1072 for (d = 0; d < coordDim; ++d) grad[fc*coordDim+d] += interpolant[fc*dim+d]*wt*detJ[q]; 1073 } 1074 vol += quadWeights[q*qNc]*detJ[q]; 1075 } 1076 fieldOffset += Nb; 1077 qc += Nc; 1078 } 1079 ierr = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr); 1080 for (fc = 0; fc < numComponents; ++fc) { 1081 for (d = 0; d < coordDim; ++d) { 1082 gradsum[fc*coordDim+d] += grad[fc*coordDim+d]; 1083 } 1084 } 1085 volsum += vol; 1086 if (debug) { 1087 ierr = PetscPrintf(PETSC_COMM_SELF, "Cell %D gradient: [", cell);CHKERRQ(ierr); 1088 for (fc = 0; fc < numComponents; ++fc) { 1089 for (d = 0; d < coordDim; ++d) { 1090 if (fc || d > 0) {ierr = PetscPrintf(PETSC_COMM_SELF, ", ");CHKERRQ(ierr);} 1091 ierr = PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc*coordDim+d]));CHKERRQ(ierr); 1092 } 1093 } 1094 ierr = PetscPrintf(PETSC_COMM_SELF, "]\n");CHKERRQ(ierr); 1095 } 1096 } 1097 for (fc = 0; fc < numComponents; ++fc) { 1098 for (d = 0; d < coordDim; ++d) gradsum[fc*coordDim+d] /= volsum; 1099 } 1100 ierr = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1101 ierr = DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES);CHKERRQ(ierr); 1102 } 1103 ierr = PetscFree6(gradsum,interpolant,coords,detJ,J,invJ);CHKERRQ(ierr); 1104 PetscFunctionReturn(0); 1105 } 1106 1107 static PetscErrorCode DMPlexComputeIntegral_Internal(DM dm, Vec X, PetscInt cStart, PetscInt cEnd, PetscScalar *cintegral, void *user) 1108 { 1109 DM dmAux = NULL; 1110 PetscDS prob, probAux; 1111 PetscSection section, sectionAux; 1112 Vec locX, locA; 1113 PetscInt dim, numCells = cEnd - cStart, cell, c, f; 1114 PetscBool useFEM = PETSC_FALSE, useFVM = PETSC_FALSE; 1115 /* DS */ 1116 PetscInt Nf, totDim, *uOff, *uOff_x, numConstants; 1117 PetscInt NfAux, totDimAux, *aOff; 1118 PetscScalar *u, *a; 1119 const PetscScalar *constants; 1120 /* Geometry */ 1121 PetscFECellGeom *cgeomFEM; 1122 DM dmGrad; 1123 Vec cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL; 1124 PetscFVCellGeom *cgeomFVM; 1125 const PetscScalar *lgrad; 1126 PetscErrorCode ierr; 1127 1128 PetscFunctionBegin; 1129 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1130 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1131 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1132 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1133 /* Determine which discretizations we have */ 1134 for (f = 0; f < Nf; ++f) { 1135 PetscObject obj; 1136 PetscClassId id; 1137 1138 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1139 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1140 if (id == PETSCFE_CLASSID) useFEM = PETSC_TRUE; 1141 if (id == PETSCFV_CLASSID) useFVM = PETSC_TRUE; 1142 } 1143 /* Get local solution with boundary values */ 1144 ierr = DMGetLocalVector(dm, &locX);CHKERRQ(ierr); 1145 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 1146 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);CHKERRQ(ierr); 1147 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);CHKERRQ(ierr); 1148 /* Read DS information */ 1149 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1150 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 1151 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 1152 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 1153 /* Read Auxiliary DS information */ 1154 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1155 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);CHKERRQ(ierr); 1156 if (dmAux) { 1157 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1158 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 1159 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1160 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1161 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 1162 } 1163 /* Allocate data arrays */ 1164 ierr = PetscCalloc1(numCells*totDim, &u);CHKERRQ(ierr); 1165 if (useFEM) {ierr = PetscMalloc1(numCells, &cgeomFEM);CHKERRQ(ierr);} 1166 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1167 /* Read out geometry */ 1168 if (useFEM) { 1169 for (cell = cStart; cell < cEnd; ++cell) { 1170 const PetscInt c = cell - cStart; 1171 1172 ierr = DMPlexComputeCellGeometryFEM(dm, cell, NULL, cgeomFEM[c].v0, cgeomFEM[c].J, cgeomFEM[c].invJ, &cgeomFEM[c].detJ);CHKERRQ(ierr); 1173 if (cgeomFEM[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeomFEM[c].detJ, c); 1174 } 1175 } 1176 if (useFVM) { 1177 PetscFV fv = NULL; 1178 Vec grad; 1179 PetscInt fStart, fEnd; 1180 PetscBool compGrad; 1181 1182 for (f = 0; f < Nf; ++f) { 1183 PetscObject obj; 1184 PetscClassId id; 1185 1186 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1187 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1188 if (id == PETSCFV_CLASSID) {fv = (PetscFV) obj; break;} 1189 } 1190 ierr = PetscFVGetComputeGradients(fv, &compGrad);CHKERRQ(ierr); 1191 ierr = PetscFVSetComputeGradients(fv, PETSC_TRUE);CHKERRQ(ierr); 1192 ierr = DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);CHKERRQ(ierr); 1193 ierr = DMPlexComputeGradientFVM(dm, fv, faceGeometryFVM, cellGeometryFVM, &dmGrad);CHKERRQ(ierr); 1194 ierr = PetscFVSetComputeGradients(fv, compGrad);CHKERRQ(ierr); 1195 ierr = VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1196 /* Reconstruct and limit cell gradients */ 1197 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1198 ierr = DMGetGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1199 ierr = DMPlexReconstructGradients_Internal(dm, fv, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);CHKERRQ(ierr); 1200 /* Communicate gradient values */ 1201 ierr = DMGetLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1202 ierr = DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1203 ierr = DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1204 ierr = DMRestoreGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1205 /* Handle non-essential (e.g. outflow) boundary values */ 1206 ierr = DMPlexInsertBoundaryValues(dm, PETSC_FALSE, locX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);CHKERRQ(ierr); 1207 ierr = VecGetArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1208 } 1209 /* Read out data from inputs */ 1210 for (c = cStart; c < cEnd; ++c) { 1211 PetscScalar *x = NULL; 1212 PetscInt i; 1213 1214 ierr = DMPlexVecGetClosure(dm, section, locX, c, NULL, &x);CHKERRQ(ierr); 1215 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 1216 ierr = DMPlexVecRestoreClosure(dm, section, locX, c, NULL, &x);CHKERRQ(ierr); 1217 if (dmAux) { 1218 ierr = DMPlexVecGetClosure(dmAux, sectionAux, locA, c, NULL, &x);CHKERRQ(ierr); 1219 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 1220 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, locA, c, NULL, &x);CHKERRQ(ierr); 1221 } 1222 } 1223 /* Do integration for each field */ 1224 for (f = 0; f < Nf; ++f) { 1225 PetscObject obj; 1226 PetscClassId id; 1227 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 1228 1229 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1230 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1231 if (id == PETSCFE_CLASSID) { 1232 PetscFE fe = (PetscFE) obj; 1233 PetscQuadrature q; 1234 PetscInt Nq, Nb; 1235 1236 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1237 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1238 ierr = PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 1239 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1240 blockSize = Nb*Nq; 1241 batchSize = numBlocks * blockSize; 1242 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1243 numChunks = numCells / (numBatches*batchSize); 1244 Ne = numChunks*numBatches*batchSize; 1245 Nr = numCells % (numBatches*batchSize); 1246 offset = numCells - Nr; 1247 ierr = PetscFEIntegrate(fe, prob, f, Ne, cgeomFEM, u, probAux, a, cintegral);CHKERRQ(ierr); 1248 ierr = PetscFEIntegrate(fe, prob, f, Nr, &cgeomFEM[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], &cintegral[offset*Nf]);CHKERRQ(ierr); 1249 } else if (id == PETSCFV_CLASSID) { 1250 PetscInt foff; 1251 PetscPointFunc obj_func; 1252 PetscScalar lint; 1253 1254 ierr = PetscDSGetObjective(prob, f, &obj_func);CHKERRQ(ierr); 1255 ierr = PetscDSGetFieldOffset(prob, f, &foff);CHKERRQ(ierr); 1256 if (obj_func) { 1257 for (c = 0; c < numCells; ++c) { 1258 PetscScalar *u_x; 1259 1260 ierr = DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);CHKERRQ(ierr); 1261 obj_func(dim, Nf, NfAux, uOff, uOff_x, &u[totDim*c+foff], NULL, u_x, aOff, NULL, &a[totDimAux*c], NULL, NULL, 0.0, cgeomFVM[c].centroid, numConstants, constants, &lint); 1262 cintegral[c*Nf+f] += PetscRealPart(lint)*cgeomFVM[c].volume; 1263 } 1264 } 1265 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 1266 } 1267 /* Cleanup data arrays */ 1268 if (useFVM) { 1269 ierr = VecRestoreArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1270 ierr = VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1271 ierr = DMRestoreLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1272 ierr = VecDestroy(&faceGeometryFVM);CHKERRQ(ierr); 1273 ierr = VecDestroy(&cellGeometryFVM);CHKERRQ(ierr); 1274 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr); 1275 } 1276 if (useFEM) {ierr = PetscFree(cgeomFEM);CHKERRQ(ierr);} 1277 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1278 ierr = PetscFree(u);CHKERRQ(ierr); 1279 /* Cleanup */ 1280 ierr = DMRestoreLocalVector(dm, &locX);CHKERRQ(ierr); 1281 PetscFunctionReturn(0); 1282 } 1283 1284 /*@ 1285 DMPlexComputeIntegralFEM - Form the integral over the domain from the global input X using pointwise functions specified by the user 1286 1287 Input Parameters: 1288 + dm - The mesh 1289 . X - Global input vector 1290 - user - The user context 1291 1292 Output Parameter: 1293 . integral - Integral for each field 1294 1295 Level: developer 1296 1297 .seealso: DMPlexComputeResidualFEM() 1298 @*/ 1299 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user) 1300 { 1301 DM_Plex *mesh = (DM_Plex *) dm->data; 1302 PetscScalar *cintegral; 1303 PetscReal *lintegral; 1304 PetscInt Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell; 1305 PetscErrorCode ierr; 1306 1307 PetscFunctionBegin; 1308 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1309 PetscValidHeaderSpecific(X, VEC_CLASSID, 2); 1310 PetscValidPointer(integral, 3); 1311 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1312 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 1313 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr); 1314 ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr); 1315 ierr = DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);CHKERRQ(ierr); 1316 cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight]; 1317 /* TODO Introduce a loop over large chunks (right now this is a single chunk) */ 1318 ierr = PetscCalloc2(Nf, &lintegral, (cEnd-cStart)*Nf, &cintegral);CHKERRQ(ierr); 1319 ierr = DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);CHKERRQ(ierr); 1320 /* Sum up values */ 1321 for (cell = cStart; cell < cEnd; ++cell) { 1322 const PetscInt c = cell - cStart; 1323 1324 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);CHKERRQ(ierr);} 1325 for (f = 0; f < Nf; ++f) lintegral[f] += PetscRealPart(cintegral[c*Nf+f]); 1326 } 1327 ierr = MPIU_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr); 1328 if (mesh->printFEM) { 1329 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Integral:");CHKERRQ(ierr); 1330 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", integral[f]);CHKERRQ(ierr);} 1331 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 1332 } 1333 ierr = PetscFree2(lintegral, cintegral);CHKERRQ(ierr); 1334 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1335 PetscFunctionReturn(0); 1336 } 1337 1338 /*@ 1339 DMPlexComputeCellwiseIntegralFEM - Form the vector of cellwise integrals F from the global input X using pointwise functions specified by the user 1340 1341 Input Parameters: 1342 + dm - The mesh 1343 . X - Global input vector 1344 - user - The user context 1345 1346 Output Parameter: 1347 . integral - Cellwise integrals for each field 1348 1349 Level: developer 1350 1351 .seealso: DMPlexComputeResidualFEM() 1352 @*/ 1353 PetscErrorCode DMPlexComputeCellwiseIntegralFEM(DM dm, Vec X, Vec F, void *user) 1354 { 1355 DM_Plex *mesh = (DM_Plex *) dm->data; 1356 DM dmF; 1357 PetscSection sectionF; 1358 PetscScalar *cintegral, *af; 1359 PetscInt Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell; 1360 PetscErrorCode ierr; 1361 1362 PetscFunctionBegin; 1363 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1364 PetscValidHeaderSpecific(X, VEC_CLASSID, 2); 1365 PetscValidHeaderSpecific(F, VEC_CLASSID, 3); 1366 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1367 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 1368 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr); 1369 ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr); 1370 ierr = DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);CHKERRQ(ierr); 1371 cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight]; 1372 /* TODO Introduce a loop over large chunks (right now this is a single chunk) */ 1373 ierr = PetscCalloc1((cEnd-cStart)*Nf, &cintegral);CHKERRQ(ierr); 1374 ierr = DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);CHKERRQ(ierr); 1375 /* Put values in F*/ 1376 ierr = VecGetDM(F, &dmF);CHKERRQ(ierr); 1377 ierr = DMGetDefaultSection(dmF, §ionF);CHKERRQ(ierr); 1378 ierr = VecGetArray(F, &af);CHKERRQ(ierr); 1379 for (cell = cStart; cell < cEnd; ++cell) { 1380 const PetscInt c = cell - cStart; 1381 PetscInt dof, off; 1382 1383 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);CHKERRQ(ierr);} 1384 ierr = PetscSectionGetDof(sectionF, cell, &dof);CHKERRQ(ierr); 1385 ierr = PetscSectionGetOffset(sectionF, cell, &off);CHKERRQ(ierr); 1386 if (dof != Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cell dofs %D != %D", dof, Nf); 1387 for (f = 0; f < Nf; ++f) af[off+f] = cintegral[c*Nf+f]; 1388 } 1389 ierr = VecRestoreArray(F, &af);CHKERRQ(ierr); 1390 ierr = PetscFree(cintegral);CHKERRQ(ierr); 1391 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1392 PetscFunctionReturn(0); 1393 } 1394 1395 /*@ 1396 DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1397 1398 Input Parameters: 1399 + dmf - The fine mesh 1400 . dmc - The coarse mesh 1401 - user - The user context 1402 1403 Output Parameter: 1404 . In - The interpolation matrix 1405 1406 Level: developer 1407 1408 .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1409 @*/ 1410 PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user) 1411 { 1412 DM_Plex *mesh = (DM_Plex *) dmc->data; 1413 const char *name = "Interpolator"; 1414 PetscDS prob; 1415 PetscFE *feRef; 1416 PetscFV *fvRef; 1417 PetscSection fsection, fglobalSection; 1418 PetscSection csection, cglobalSection; 1419 PetscScalar *elemMat; 1420 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c; 1421 PetscInt cTotDim, rTotDim = 0; 1422 PetscErrorCode ierr; 1423 1424 PetscFunctionBegin; 1425 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1426 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1427 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1428 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1429 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1430 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1431 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1432 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1433 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1434 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1435 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr); 1436 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1437 for (f = 0; f < Nf; ++f) { 1438 PetscObject obj; 1439 PetscClassId id; 1440 PetscInt rNb = 0, Nc = 0; 1441 1442 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1443 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1444 if (id == PETSCFE_CLASSID) { 1445 PetscFE fe = (PetscFE) obj; 1446 1447 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1448 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1449 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1450 } else if (id == PETSCFV_CLASSID) { 1451 PetscFV fv = (PetscFV) obj; 1452 PetscDualSpace Q; 1453 1454 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1455 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1456 ierr = PetscDualSpaceGetDimension(Q, &rNb);CHKERRQ(ierr); 1457 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1458 } 1459 rTotDim += rNb; 1460 } 1461 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1462 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 1463 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 1464 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1465 PetscDualSpace Qref; 1466 PetscQuadrature f; 1467 const PetscReal *qpoints, *qweights; 1468 PetscReal *points; 1469 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1470 1471 /* Compose points from all dual basis functionals */ 1472 if (feRef[fieldI]) { 1473 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1474 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 1475 } else { 1476 ierr = PetscFVGetDualSpace(fvRef[fieldI], &Qref);CHKERRQ(ierr); 1477 ierr = PetscFVGetNumComponents(fvRef[fieldI], &Nc);CHKERRQ(ierr); 1478 } 1479 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1480 for (i = 0; i < fpdim; ++i) { 1481 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1482 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1483 npoints += Np; 1484 } 1485 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1486 for (i = 0, k = 0; i < fpdim; ++i) { 1487 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1488 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1489 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1490 } 1491 1492 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1493 PetscObject obj; 1494 PetscClassId id; 1495 PetscReal *B; 1496 PetscInt NcJ = 0, cpdim = 0, j, qNc; 1497 1498 ierr = PetscDSGetDiscretization(prob, fieldJ, &obj);CHKERRQ(ierr); 1499 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1500 if (id == PETSCFE_CLASSID) { 1501 PetscFE fe = (PetscFE) obj; 1502 1503 /* Evaluate basis at points */ 1504 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 1505 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1506 /* For now, fields only interpolate themselves */ 1507 if (fieldI == fieldJ) { 1508 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); 1509 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1510 for (i = 0, k = 0; i < fpdim; ++i) { 1511 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1512 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);CHKERRQ(ierr); 1513 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); 1514 for (p = 0; p < Np; ++p, ++k) { 1515 for (j = 0; j < cpdim; ++j) { 1516 /* NOTE: This is not quite right, unless fpdim == number of fine grid functional quad points */ 1517 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p*qNc+c]; 1518 } 1519 } 1520 } 1521 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1522 } 1523 } else if (id == PETSCFV_CLASSID) { 1524 PetscFV fv = (PetscFV) obj; 1525 1526 /* Evaluate constant function at points */ 1527 ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr); 1528 cpdim = 1; 1529 /* For now, fields only interpolate themselves */ 1530 if (fieldI == fieldJ) { 1531 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); 1532 for (i = 0, k = 0; i < fpdim; ++i) { 1533 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1534 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);CHKERRQ(ierr); 1535 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); 1536 for (p = 0; p < Np; ++p, ++k) { 1537 for (j = 0; j < cpdim; ++j) { 1538 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p*qNc+c]; 1539 } 1540 } 1541 } 1542 } 1543 } 1544 offsetJ += cpdim*NcJ; 1545 } 1546 offsetI += fpdim*Nc; 1547 ierr = PetscFree(points);CHKERRQ(ierr); 1548 } 1549 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1550 /* Preallocate matrix */ 1551 { 1552 Mat preallocator; 1553 PetscScalar *vals; 1554 PetscInt *cellCIndices, *cellFIndices; 1555 PetscInt locRows, locCols, cell; 1556 1557 ierr = MatGetLocalSize(In, &locRows, &locCols);CHKERRQ(ierr); 1558 ierr = MatCreate(PetscObjectComm((PetscObject) In), &preallocator);CHKERRQ(ierr); 1559 ierr = MatSetType(preallocator, MATPREALLOCATOR);CHKERRQ(ierr); 1560 ierr = MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1561 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 1562 ierr = PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1563 for (cell = cStart; cell < cEnd; ++cell) { 1564 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1565 ierr = MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);CHKERRQ(ierr); 1566 } 1567 ierr = PetscFree3(vals,cellCIndices,cellFIndices);CHKERRQ(ierr); 1568 ierr = MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1569 ierr = MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1570 ierr = MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);CHKERRQ(ierr); 1571 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 1572 } 1573 /* Fill matrix */ 1574 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1575 for (c = cStart; c < cEnd; ++c) { 1576 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1577 } 1578 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1579 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1580 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1581 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1582 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1583 if (mesh->printFEM) { 1584 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1585 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1586 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1587 } 1588 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1589 PetscFunctionReturn(0); 1590 } 1591 1592 PetscErrorCode DMPlexComputeMassMatrixNested(DM dmc, DM dmf, Mat mass, void *user) 1593 { 1594 SETERRQ(PetscObjectComm((PetscObject) dmc), PETSC_ERR_SUP, "Laziness"); 1595 } 1596 1597 /*@ 1598 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 1599 1600 Input Parameters: 1601 + dmf - The fine mesh 1602 . dmc - The coarse mesh 1603 - user - The user context 1604 1605 Output Parameter: 1606 . In - The interpolation matrix 1607 1608 Level: developer 1609 1610 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1611 @*/ 1612 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 1613 { 1614 DM_Plex *mesh = (DM_Plex *) dmf->data; 1615 const char *name = "Interpolator"; 1616 PetscDS prob; 1617 PetscSection fsection, csection, globalFSection, globalCSection; 1618 PetscHashJK ht; 1619 PetscLayout rLayout; 1620 PetscInt *dnz, *onz; 1621 PetscInt locRows, rStart, rEnd; 1622 PetscReal *x, *v0, *J, *invJ, detJ; 1623 PetscReal *v0c, *Jc, *invJc, detJc; 1624 PetscScalar *elemMat; 1625 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1626 PetscErrorCode ierr; 1627 1628 PetscFunctionBegin; 1629 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1630 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1631 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1632 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1633 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1634 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1635 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1636 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1637 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1638 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1639 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1640 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1641 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1642 ierr = PetscMalloc1(totDim, &elemMat);CHKERRQ(ierr); 1643 1644 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1645 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1646 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1647 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1648 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1649 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1650 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1651 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1652 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1653 for (field = 0; field < Nf; ++field) { 1654 PetscObject obj; 1655 PetscClassId id; 1656 PetscDualSpace Q = NULL; 1657 PetscQuadrature f; 1658 const PetscReal *qpoints; 1659 PetscInt Nc, Np, fpdim, i, d; 1660 1661 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1662 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1663 if (id == PETSCFE_CLASSID) { 1664 PetscFE fe = (PetscFE) obj; 1665 1666 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1667 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1668 } else if (id == PETSCFV_CLASSID) { 1669 PetscFV fv = (PetscFV) obj; 1670 1671 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1672 Nc = 1; 1673 } 1674 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1675 /* For each fine grid cell */ 1676 for (cell = cStart; cell < cEnd; ++cell) { 1677 PetscInt *findices, *cindices; 1678 PetscInt numFIndices, numCIndices; 1679 1680 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1681 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1682 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim); 1683 for (i = 0; i < fpdim; ++i) { 1684 Vec pointVec; 1685 PetscScalar *pV; 1686 PetscSF coarseCellSF = NULL; 1687 const PetscSFNode *coarseCells; 1688 PetscInt numCoarseCells, q, c; 1689 1690 /* Get points from the dual basis functional quadrature */ 1691 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1692 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1693 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1694 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1695 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1696 for (q = 0; q < Np; ++q) { 1697 /* Transform point to real space */ 1698 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1699 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1700 } 1701 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1702 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1703 /* OPT: Pack all quad points from fine cell */ 1704 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1705 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 1706 /* Update preallocation info */ 1707 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1708 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1709 { 1710 PetscHashJKKey key; 1711 PetscHashJKIter missing, iter; 1712 1713 key.j = findices[i]; 1714 if (key.j >= 0) { 1715 /* Get indices for coarse elements */ 1716 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1717 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1718 for (c = 0; c < numCIndices; ++c) { 1719 key.k = cindices[c]; 1720 if (key.k < 0) continue; 1721 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1722 if (missing) { 1723 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1724 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1725 else ++onz[key.j-rStart]; 1726 } 1727 } 1728 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1729 } 1730 } 1731 } 1732 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1733 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1734 } 1735 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1736 } 1737 } 1738 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1739 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1740 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1741 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1742 for (field = 0; field < Nf; ++field) { 1743 PetscObject obj; 1744 PetscClassId id; 1745 PetscDualSpace Q = NULL; 1746 PetscQuadrature f; 1747 const PetscReal *qpoints, *qweights; 1748 PetscInt Nc, qNc, Np, fpdim, i, d; 1749 1750 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1751 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1752 if (id == PETSCFE_CLASSID) { 1753 PetscFE fe = (PetscFE) obj; 1754 1755 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1756 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1757 } else if (id == PETSCFV_CLASSID) { 1758 PetscFV fv = (PetscFV) obj; 1759 1760 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1761 Nc = 1; 1762 } 1763 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1764 /* For each fine grid cell */ 1765 for (cell = cStart; cell < cEnd; ++cell) { 1766 PetscInt *findices, *cindices; 1767 PetscInt numFIndices, numCIndices; 1768 1769 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1770 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1771 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim); 1772 for (i = 0; i < fpdim; ++i) { 1773 Vec pointVec; 1774 PetscScalar *pV; 1775 PetscSF coarseCellSF = NULL; 1776 const PetscSFNode *coarseCells; 1777 PetscInt numCoarseCells, cpdim, q, c, j; 1778 1779 /* Get points from the dual basis functional quadrature */ 1780 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1781 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);CHKERRQ(ierr); 1782 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); 1783 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1784 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1785 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1786 for (q = 0; q < Np; ++q) { 1787 /* Transform point to real space */ 1788 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1789 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1790 } 1791 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1792 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1793 /* OPT: Read this out from preallocation information */ 1794 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1795 /* Update preallocation info */ 1796 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1797 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1798 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1799 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1800 PetscReal pVReal[3]; 1801 1802 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1803 /* Transform points from real space to coarse reference space */ 1804 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 1805 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 1806 CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x); 1807 1808 if (id == PETSCFE_CLASSID) { 1809 PetscFE fe = (PetscFE) obj; 1810 PetscReal *B; 1811 1812 /* Evaluate coarse basis on contained point */ 1813 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1814 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 1815 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 1816 /* Get elemMat entries by multiplying by weight */ 1817 for (j = 0; j < cpdim; ++j) { 1818 for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*qweights[ccell*qNc + c]; 1819 } 1820 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1821 } else { 1822 cpdim = 1; 1823 for (j = 0; j < cpdim; ++j) { 1824 for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c]; 1825 } 1826 } 1827 /* Update interpolator */ 1828 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 1829 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 1830 ierr = MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1831 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1832 } 1833 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1834 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1835 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1836 } 1837 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1838 } 1839 } 1840 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 1841 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 1842 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1843 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1844 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1845 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1846 PetscFunctionReturn(0); 1847 } 1848 1849 /*@ 1850 DMPlexComputeMassMatrixGeneral - Form the local portion of the mass matrix M from the coarse DM to a non-nested fine DM. 1851 1852 Input Parameters: 1853 + dmf - The fine mesh 1854 . dmc - The coarse mesh 1855 - user - The user context 1856 1857 Output Parameter: 1858 . mass - The mass matrix 1859 1860 Level: developer 1861 1862 .seealso: DMPlexComputeMassMatrixNested(), DMPlexComputeInterpolatorNested(), DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1863 @*/ 1864 PetscErrorCode DMPlexComputeMassMatrixGeneral(DM dmc, DM dmf, Mat mass, void *user) 1865 { 1866 DM_Plex *mesh = (DM_Plex *) dmf->data; 1867 const char *name = "Mass Matrix"; 1868 PetscDS prob; 1869 PetscSection fsection, csection, globalFSection, globalCSection; 1870 PetscHashJK ht; 1871 PetscLayout rLayout; 1872 PetscInt *dnz, *onz; 1873 PetscInt locRows, rStart, rEnd; 1874 PetscReal *x, *v0, *J, *invJ, detJ; 1875 PetscReal *v0c, *Jc, *invJc, detJc; 1876 PetscScalar *elemMat; 1877 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1878 PetscErrorCode ierr; 1879 1880 PetscFunctionBegin; 1881 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1882 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1883 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1884 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1885 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1886 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1887 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1888 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1889 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1890 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1891 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1892 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1893 ierr = PetscMalloc1(totDim, &elemMat);CHKERRQ(ierr); 1894 1895 ierr = MatGetLocalSize(mass, &locRows, NULL);CHKERRQ(ierr); 1896 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) mass), &rLayout);CHKERRQ(ierr); 1897 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1898 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1899 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1900 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1901 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1902 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1903 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1904 for (field = 0; field < Nf; ++field) { 1905 PetscObject obj; 1906 PetscClassId id; 1907 PetscQuadrature quad; 1908 const PetscReal *qpoints; 1909 PetscInt Nq, Nc, i, d; 1910 1911 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1912 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1913 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetQuadrature((PetscFE) obj, &quad);CHKERRQ(ierr);} 1914 else {ierr = PetscFVGetQuadrature((PetscFV) obj, &quad);CHKERRQ(ierr);} 1915 ierr = PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, NULL);CHKERRQ(ierr); 1916 /* For each fine grid cell */ 1917 for (cell = cStart; cell < cEnd; ++cell) { 1918 Vec pointVec; 1919 PetscScalar *pV; 1920 PetscSF coarseCellSF = NULL; 1921 const PetscSFNode *coarseCells; 1922 PetscInt numCoarseCells, q, c; 1923 PetscInt *findices, *cindices; 1924 PetscInt numFIndices, numCIndices; 1925 1926 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1927 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1928 /* Get points from the quadrature */ 1929 ierr = VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);CHKERRQ(ierr); 1930 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1931 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1932 for (q = 0; q < Nq; ++q) { 1933 /* Transform point to real space */ 1934 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1935 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1936 } 1937 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1938 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1939 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1940 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 1941 /* Update preallocation info */ 1942 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1943 if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1944 { 1945 PetscHashJKKey key; 1946 PetscHashJKIter missing, iter; 1947 1948 for (i = 0; i < numFIndices; ++i) { 1949 key.j = findices[i]; 1950 if (key.j >= 0) { 1951 /* Get indices for coarse elements */ 1952 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1953 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1954 for (c = 0; c < numCIndices; ++c) { 1955 key.k = cindices[c]; 1956 if (key.k < 0) continue; 1957 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1958 if (missing) { 1959 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1960 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1961 else ++onz[key.j-rStart]; 1962 } 1963 } 1964 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1965 } 1966 } 1967 } 1968 } 1969 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1970 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1971 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1972 } 1973 } 1974 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1975 ierr = MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1976 ierr = MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1977 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1978 for (field = 0; field < Nf; ++field) { 1979 PetscObject obj; 1980 PetscClassId id; 1981 PetscQuadrature quad; 1982 PetscReal *Bfine; 1983 const PetscReal *qpoints, *qweights; 1984 PetscInt Nq, Nc, i, d; 1985 1986 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1987 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1988 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetQuadrature((PetscFE) obj, &quad);CHKERRQ(ierr);ierr = PetscFEGetDefaultTabulation((PetscFE) obj, &Bfine, NULL, NULL);CHKERRQ(ierr);} 1989 else {ierr = PetscFVGetQuadrature((PetscFV) obj, &quad);CHKERRQ(ierr);} 1990 ierr = PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, &qweights);CHKERRQ(ierr); 1991 /* For each fine grid cell */ 1992 for (cell = cStart; cell < cEnd; ++cell) { 1993 Vec pointVec; 1994 PetscScalar *pV; 1995 PetscSF coarseCellSF = NULL; 1996 const PetscSFNode *coarseCells; 1997 PetscInt numCoarseCells, cpdim, q, c, j; 1998 PetscInt *findices, *cindices; 1999 PetscInt numFIndices, numCIndices; 2000 2001 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2002 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 2003 /* Get points from the quadrature */ 2004 ierr = VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);CHKERRQ(ierr); 2005 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 2006 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2007 for (q = 0; q < Nq; ++q) { 2008 /* Transform point to real space */ 2009 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 2010 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 2011 } 2012 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2013 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 2014 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 2015 /* Update matrix */ 2016 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 2017 if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 2018 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2019 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 2020 PetscReal pVReal[3]; 2021 2022 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2023 /* Transform points from real space to coarse reference space */ 2024 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 2025 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 2026 CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x); 2027 2028 if (id == PETSCFE_CLASSID) { 2029 PetscFE fe = (PetscFE) obj; 2030 PetscReal *B; 2031 2032 /* Evaluate coarse basis on contained point */ 2033 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 2034 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 2035 /* Get elemMat entries by multiplying by weight */ 2036 for (i = 0; i < numFIndices; ++i) { 2037 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 2038 for (j = 0; j < cpdim; ++j) { 2039 for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*Bfine[(ccell*numFIndices + i)*Nc + c]*qweights[ccell*Nc + c]*detJ; 2040 } 2041 /* Update interpolator */ 2042 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 2043 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 2044 ierr = MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);CHKERRQ(ierr); 2045 } 2046 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 2047 } else { 2048 cpdim = 1; 2049 for (i = 0; i < numFIndices; ++i) { 2050 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 2051 for (j = 0; j < cpdim; ++j) { 2052 for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*1.0*qweights[ccell*Nc + c]*detJ; 2053 } 2054 /* Update interpolator */ 2055 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 2056 ierr = PetscPrintf(PETSC_COMM_SELF, "Nq: %d %d Nf: %d %d Nc: %d %d\n", ccell, Nq, i, numFIndices, j, numCIndices);CHKERRQ(ierr); 2057 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 2058 ierr = MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);CHKERRQ(ierr); 2059 } 2060 } 2061 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2062 } 2063 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2064 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 2065 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 2066 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2067 } 2068 } 2069 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 2070 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 2071 ierr = PetscFree(elemMat);CHKERRQ(ierr); 2072 ierr = MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2073 ierr = MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2074 PetscFunctionReturn(0); 2075 } 2076 2077 /*@ 2078 DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns 2079 2080 Input Parameters: 2081 + dmc - The coarse mesh 2082 - dmf - The fine mesh 2083 - user - The user context 2084 2085 Output Parameter: 2086 . sc - The mapping 2087 2088 Level: developer 2089 2090 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 2091 @*/ 2092 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 2093 { 2094 PetscDS prob; 2095 PetscFE *feRef; 2096 PetscFV *fvRef; 2097 Vec fv, cv; 2098 IS fis, cis; 2099 PetscSection fsection, fglobalSection, csection, cglobalSection; 2100 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 2101 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m; 2102 PetscErrorCode ierr; 2103 2104 PetscFunctionBegin; 2105 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2106 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 2107 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 2108 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 2109 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 2110 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 2111 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 2112 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 2113 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 2114 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 2115 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 2116 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 2117 for (f = 0; f < Nf; ++f) { 2118 PetscObject obj; 2119 PetscClassId id; 2120 PetscInt fNb = 0, Nc = 0; 2121 2122 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 2123 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2124 if (id == PETSCFE_CLASSID) { 2125 PetscFE fe = (PetscFE) obj; 2126 2127 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 2128 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 2129 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 2130 } else if (id == PETSCFV_CLASSID) { 2131 PetscFV fv = (PetscFV) obj; 2132 PetscDualSpace Q; 2133 2134 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 2135 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 2136 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 2137 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 2138 } 2139 fTotDim += fNb*Nc; 2140 } 2141 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 2142 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 2143 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 2144 PetscFE feC; 2145 PetscFV fvC; 2146 PetscDualSpace QF, QC; 2147 PetscInt NcF, NcC, fpdim, cpdim; 2148 2149 if (feRef[field]) { 2150 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 2151 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 2152 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 2153 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 2154 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 2155 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 2156 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 2157 } else { 2158 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 2159 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 2160 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 2161 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 2162 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 2163 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 2164 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 2165 } 2166 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); 2167 for (c = 0; c < cpdim; ++c) { 2168 PetscQuadrature cfunc; 2169 const PetscReal *cqpoints; 2170 PetscInt NpC; 2171 PetscBool found = PETSC_FALSE; 2172 2173 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 2174 ierr = PetscQuadratureGetData(cfunc, NULL, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 2175 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 2176 for (f = 0; f < fpdim; ++f) { 2177 PetscQuadrature ffunc; 2178 const PetscReal *fqpoints; 2179 PetscReal sum = 0.0; 2180 PetscInt NpF, comp; 2181 2182 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 2183 ierr = PetscQuadratureGetData(ffunc, NULL, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 2184 if (NpC != NpF) continue; 2185 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 2186 if (sum > 1.0e-9) continue; 2187 for (comp = 0; comp < NcC; ++comp) { 2188 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 2189 } 2190 found = PETSC_TRUE; 2191 break; 2192 } 2193 if (!found) { 2194 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 2195 if (fvRef[field]) { 2196 PetscInt comp; 2197 for (comp = 0; comp < NcC; ++comp) { 2198 cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp; 2199 } 2200 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 2201 } 2202 } 2203 offsetC += cpdim*NcC; 2204 offsetF += fpdim*NcF; 2205 } 2206 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 2207 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 2208 2209 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 2210 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 2211 ierr = VecGetOwnershipRange(cv, &startC, &endC);CHKERRQ(ierr); 2212 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 2213 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 2214 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 2215 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 2216 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 2217 for (c = cStart; c < cEnd; ++c) { 2218 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 2219 for (d = 0; d < cTotDim; ++d) { 2220 if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue; 2221 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]]); 2222 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 2223 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 2224 } 2225 } 2226 ierr = PetscFree(cmap);CHKERRQ(ierr); 2227 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 2228 2229 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 2230 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 2231 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 2232 ierr = ISDestroy(&cis);CHKERRQ(ierr); 2233 ierr = ISDestroy(&fis);CHKERRQ(ierr); 2234 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 2235 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 2236 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2237 PetscFunctionReturn(0); 2238 } 2239