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