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)(void); 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, 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, &Nq, NULL, &quadWeights);CHKERRQ(ierr); 511 ierr = PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);CHKERRQ(ierr); 512 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 513 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 514 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 515 for (c = cStart; c < cEnd; ++c) { 516 PetscScalar *x = NULL; 517 PetscReal elemDiff = 0.0; 518 519 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);CHKERRQ(ierr); 520 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 521 522 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 523 PetscObject obj; 524 PetscClassId id; 525 void * const ctx = ctxs ? ctxs[field] : NULL; 526 PetscInt Nb, Nc, q, fc; 527 528 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 529 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 530 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 531 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 532 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 533 if (debug) { 534 char title[1024]; 535 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 536 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 537 } 538 for (q = 0; q < Nq; ++q) { 539 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); 540 ierr = (*funcs[field])(coordDim, time, &coords[coordDim * q], Nc, funcVal, ctx); 541 if (ierr) { 542 PetscErrorCode ierr2; 543 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 544 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 545 ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2); 546 CHKERRQ(ierr); 547 } 548 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 549 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 550 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 551 for (fc = 0; fc < Nc; ++fc) { 552 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ[q]);CHKERRQ(ierr);} 553 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ[q]; 554 } 555 } 556 fieldOffset += Nb*Nc; 557 } 558 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 559 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 560 localDiff += elemDiff; 561 } 562 ierr = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr); 563 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 564 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 565 *diff = PetscSqrtReal(*diff); 566 PetscFunctionReturn(0); 567 } 568 569 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) 570 { 571 const PetscInt debug = 0; 572 PetscSection section; 573 PetscQuadrature quad; 574 Vec localX; 575 PetscScalar *funcVal, *interpolantVec; 576 PetscReal *coords, *realSpaceDer, *J, *invJ, *detJ; 577 PetscReal localDiff = 0.0; 578 PetscInt dim, coordDim, Nq = 0, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset, comp; 579 PetscErrorCode ierr; 580 581 PetscFunctionBegin; 582 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 583 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 584 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 585 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 586 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 587 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 588 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 589 for (field = 0; field < numFields; ++field) { 590 PetscFE fe; 591 PetscInt Nc; 592 593 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 594 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 595 ierr = PetscQuadratureGetData(quad, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 596 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 597 numComponents += Nc; 598 } 599 /* ierr = DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 600 ierr = PetscMalloc7(numComponents,&funcVal,coordDim*Nq,&coords,coordDim*Nq,&realSpaceDer,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ,coordDim,&interpolantVec,Nq,&detJ);CHKERRQ(ierr); 601 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 602 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 603 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 604 for (c = cStart; c < cEnd; ++c) { 605 PetscScalar *x = NULL; 606 PetscReal elemDiff = 0.0; 607 608 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, invJ, detJ);CHKERRQ(ierr); 609 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 610 611 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 612 PetscFE fe; 613 void * const ctx = ctxs ? ctxs[field] : NULL; 614 const PetscReal *quadPoints, *quadWeights; 615 PetscReal *basisDer; 616 PetscInt numQuadPoints, Nb, Ncomp, q, d, fc, f, g; 617 618 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 619 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 620 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 621 ierr = PetscFEGetNumComponents(fe, &Ncomp);CHKERRQ(ierr); 622 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 623 if (debug) { 624 char title[1024]; 625 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 626 ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr); 627 } 628 for (q = 0; q < numQuadPoints; ++q) { 629 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); 630 ierr = (*funcs[field])(coordDim, time, &coords[q*coordDim], n, numFields, funcVal, ctx); 631 if (ierr) { 632 PetscErrorCode ierr2; 633 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 634 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 635 ierr2 = PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolantVec,detJ);CHKERRQ(ierr2); 636 CHKERRQ(ierr); 637 } 638 for (fc = 0; fc < Ncomp; ++fc) { 639 PetscScalar interpolant = 0.0; 640 641 for (d = 0; d < coordDim; ++d) interpolantVec[d] = 0.0; 642 for (f = 0; f < Nb; ++f) { 643 const PetscInt fidx = f*Ncomp+fc; 644 645 for (d = 0; d < coordDim; ++d) { 646 realSpaceDer[d] = 0.0; 647 for (g = 0; g < dim; ++g) { 648 realSpaceDer[d] += invJ[coordDim * coordDim * q + g*coordDim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g]; 649 } 650 interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d]; 651 } 652 } 653 for (d = 0; d < coordDim; ++d) interpolant += interpolantVec[d]*n[d]; 654 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ[q]);CHKERRQ(ierr);} 655 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ[q]; 656 } 657 } 658 comp += Ncomp; 659 fieldOffset += Nb*Ncomp; 660 } 661 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 662 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 663 localDiff += elemDiff; 664 } 665 ierr = PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolantVec,detJ);CHKERRQ(ierr); 666 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 667 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 668 *diff = PetscSqrtReal(*diff); 669 PetscFunctionReturn(0); 670 } 671 672 PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 673 { 674 const PetscInt debug = 0; 675 PetscSection section; 676 PetscQuadrature quad; 677 Vec localX; 678 PetscScalar *funcVal, *interpolant; 679 PetscReal *coords, *detJ, *J; 680 PetscReal *localDiff; 681 const PetscReal *quadPoints, *quadWeights; 682 PetscInt dim, coordDim, numFields, numComponents = 0, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset; 683 PetscErrorCode ierr; 684 685 PetscFunctionBegin; 686 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 687 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 688 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 689 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 690 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 691 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 692 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 693 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 694 for (field = 0; field < numFields; ++field) { 695 PetscObject obj; 696 PetscClassId id; 697 PetscInt Nc; 698 699 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 700 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 701 if (id == PETSCFE_CLASSID) { 702 PetscFE fe = (PetscFE) obj; 703 704 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 705 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 706 } else if (id == PETSCFV_CLASSID) { 707 PetscFV fv = (PetscFV) obj; 708 709 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 710 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 711 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 712 numComponents += Nc; 713 } 714 ierr = PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 715 ierr = PetscCalloc6(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);CHKERRQ(ierr); 716 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 717 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 718 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 719 for (c = cStart; c < cEnd; ++c) { 720 PetscScalar *x = NULL; 721 722 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);CHKERRQ(ierr); 723 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 724 725 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 726 PetscObject obj; 727 PetscClassId id; 728 void * const ctx = ctxs ? ctxs[field] : NULL; 729 PetscInt Nb, Nc, q, fc; 730 731 PetscReal elemDiff = 0.0; 732 733 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 734 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 735 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 736 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 737 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 738 if (debug) { 739 char title[1024]; 740 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 741 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 742 } 743 for (q = 0; q < Nq; ++q) { 744 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); 745 ierr = (*funcs[field])(coordDim, time, &coords[coordDim*q], numFields, funcVal, ctx); 746 if (ierr) { 747 PetscErrorCode ierr2; 748 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 749 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 750 ierr2 = PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2); 751 CHKERRQ(ierr); 752 } 753 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 754 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 755 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 756 for (fc = 0; fc < Nc; ++fc) { 757 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]))*quadWeights[q]*detJ[q]);CHKERRQ(ierr);} 758 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ[q]; 759 } 760 } 761 fieldOffset += Nb*Nc; 762 localDiff[field] += elemDiff; 763 } 764 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 765 } 766 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 767 ierr = MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 768 for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]); 769 ierr = PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr); 770 PetscFunctionReturn(0); 771 } 772 773 /*@C 774 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. 775 776 Input Parameters: 777 + dm - The DM 778 . time - The time 779 . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation 780 . ctxs - Optional array of contexts to pass to each function, or NULL. 781 - X - The coefficient vector u_h 782 783 Output Parameter: 784 . D - A Vec which holds the difference ||u - u_h||_2 for each cell 785 786 Level: developer 787 788 .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff() 789 @*/ 790 PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D) 791 { 792 PetscSection section; 793 PetscQuadrature quad; 794 Vec localX; 795 PetscScalar *funcVal, *interpolant; 796 PetscReal *coords, *detJ, *J; 797 const PetscReal *quadPoints, *quadWeights; 798 PetscInt dim, coordDim, numFields, numComponents = 0, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset; 799 PetscErrorCode ierr; 800 801 PetscFunctionBegin; 802 ierr = VecSet(D, 0.0);CHKERRQ(ierr); 803 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 804 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 805 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 806 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 807 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 808 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 809 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 810 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 811 for (field = 0; field < numFields; ++field) { 812 PetscObject obj; 813 PetscClassId id; 814 PetscInt Nc; 815 816 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 817 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 818 if (id == PETSCFE_CLASSID) { 819 PetscFE fe = (PetscFE) obj; 820 821 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 822 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 823 } else if (id == PETSCFV_CLASSID) { 824 PetscFV fv = (PetscFV) obj; 825 826 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 827 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 828 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 829 numComponents += Nc; 830 } 831 ierr = PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 832 ierr = PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);CHKERRQ(ierr); 833 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 834 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 835 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 836 for (c = cStart; c < cEnd; ++c) { 837 PetscScalar *x = NULL; 838 PetscScalar elemDiff = 0.0; 839 840 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);CHKERRQ(ierr); 841 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 842 843 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 844 PetscObject obj; 845 PetscClassId id; 846 void * const ctx = ctxs ? ctxs[field] : NULL; 847 PetscInt Nb, Nc, q, fc; 848 849 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 850 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 851 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 852 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 853 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 854 if (funcs[field]) { 855 for (q = 0; q < Nq; ++q) { 856 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); 857 ierr = (*funcs[field])(dim, time, &coords[q*coordDim], Nc, funcVal, ctx); 858 if (ierr) { 859 PetscErrorCode ierr2; 860 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 861 ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr); 862 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 863 CHKERRQ(ierr); 864 } 865 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 866 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 867 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 868 for (fc = 0; fc < Nc; ++fc) { 869 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ[q]; 870 } 871 } 872 } 873 fieldOffset += Nb*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, &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*Nc; 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, &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, &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; 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, &Np, NULL, &qweights);CHKERRQ(ierr); 1181 for (p = 0; p < Np; ++p, ++k) { 1182 for (j = 0; j < cpdim; ++j) { 1183 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p]; 1184 } 1185 } 1186 } 1187 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1188 } 1189 } else if (id == PETSCFV_CLASSID) { 1190 PetscFV fv = (PetscFV) obj; 1191 1192 /* Evaluate constant function at points */ 1193 ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr); 1194 cpdim = 1; 1195 /* For now, fields only interpolate themselves */ 1196 if (fieldI == fieldJ) { 1197 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); 1198 for (i = 0, k = 0; i < fpdim; ++i) { 1199 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1200 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1201 for (p = 0; p < Np; ++p, ++k) { 1202 for (j = 0; j < cpdim; ++j) { 1203 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p]; 1204 } 1205 } 1206 } 1207 } 1208 } 1209 offsetJ += cpdim*NcJ; 1210 } 1211 offsetI += fpdim*Nc; 1212 ierr = PetscFree(points);CHKERRQ(ierr); 1213 } 1214 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1215 /* Preallocate matrix */ 1216 { 1217 Mat preallocator; 1218 PetscScalar *vals; 1219 PetscInt *cellCIndices, *cellFIndices; 1220 PetscInt locRows, locCols, cell; 1221 1222 ierr = MatGetLocalSize(In, &locRows, &locCols);CHKERRQ(ierr); 1223 ierr = MatCreate(PetscObjectComm((PetscObject) In), &preallocator);CHKERRQ(ierr); 1224 ierr = MatSetType(preallocator, MATPREALLOCATOR);CHKERRQ(ierr); 1225 ierr = MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1226 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 1227 ierr = PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1228 for (cell = cStart; cell < cEnd; ++cell) { 1229 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1230 ierr = MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);CHKERRQ(ierr); 1231 } 1232 ierr = PetscFree3(vals,cellCIndices,cellFIndices);CHKERRQ(ierr); 1233 ierr = MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1234 ierr = MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1235 ierr = MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);CHKERRQ(ierr); 1236 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 1237 } 1238 /* Fill matrix */ 1239 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1240 for (c = cStart; c < cEnd; ++c) { 1241 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1242 } 1243 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1244 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1245 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1246 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1247 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1248 if (mesh->printFEM) { 1249 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1250 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1251 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1252 } 1253 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1254 PetscFunctionReturn(0); 1255 } 1256 1257 /*@ 1258 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 1259 1260 Input Parameters: 1261 + dmf - The fine mesh 1262 . dmc - The coarse mesh 1263 - user - The user context 1264 1265 Output Parameter: 1266 . In - The interpolation matrix 1267 1268 Level: developer 1269 1270 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1271 @*/ 1272 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 1273 { 1274 DM_Plex *mesh = (DM_Plex *) dmf->data; 1275 const char *name = "Interpolator"; 1276 PetscDS prob; 1277 PetscSection fsection, csection, globalFSection, globalCSection; 1278 PetscHashJK ht; 1279 PetscLayout rLayout; 1280 PetscInt *dnz, *onz; 1281 PetscInt locRows, rStart, rEnd; 1282 PetscReal *x, *v0, *J, *invJ, detJ; 1283 PetscReal *v0c, *Jc, *invJc, detJc; 1284 PetscScalar *elemMat; 1285 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1286 PetscErrorCode ierr; 1287 1288 PetscFunctionBegin; 1289 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1290 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1291 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1292 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1293 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1294 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1295 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1296 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1297 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1298 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1299 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1300 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1301 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1302 ierr = PetscMalloc1(totDim*totDim, &elemMat);CHKERRQ(ierr); 1303 1304 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1305 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1306 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1307 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1308 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1309 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1310 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1311 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1312 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1313 for (field = 0; field < Nf; ++field) { 1314 PetscObject obj; 1315 PetscClassId id; 1316 PetscDualSpace Q = NULL; 1317 PetscQuadrature f; 1318 const PetscReal *qpoints; 1319 PetscInt Nc, Np, fpdim, i, d; 1320 1321 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1322 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1323 if (id == PETSCFE_CLASSID) { 1324 PetscFE fe = (PetscFE) obj; 1325 1326 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1327 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1328 } else if (id == PETSCFV_CLASSID) { 1329 PetscFV fv = (PetscFV) obj; 1330 1331 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1332 Nc = 1; 1333 } 1334 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1335 /* For each fine grid cell */ 1336 for (cell = cStart; cell < cEnd; ++cell) { 1337 PetscInt *findices, *cindices; 1338 PetscInt numFIndices, numCIndices; 1339 1340 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1341 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1342 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1343 for (i = 0; i < fpdim; ++i) { 1344 Vec pointVec; 1345 PetscScalar *pV; 1346 PetscSF coarseCellSF = NULL; 1347 const PetscSFNode *coarseCells; 1348 PetscInt numCoarseCells, q, r, c; 1349 1350 /* Get points from the dual basis functional quadrature */ 1351 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1352 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1353 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1354 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1355 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1356 for (q = 0; q < Np; ++q) { 1357 /* Transform point to real space */ 1358 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1359 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1360 } 1361 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1362 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1363 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1364 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 1365 /* Update preallocation info */ 1366 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1367 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1368 for (r = 0; r < Nc; ++r) { 1369 PetscHashJKKey key; 1370 PetscHashJKIter missing, iter; 1371 1372 key.j = findices[i*Nc+r]; 1373 if (key.j < 0) continue; 1374 /* Get indices for coarse elements */ 1375 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1376 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1377 for (c = 0; c < numCIndices; ++c) { 1378 key.k = cindices[c]; 1379 if (key.k < 0) continue; 1380 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1381 if (missing) { 1382 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1383 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1384 else ++onz[key.j-rStart]; 1385 } 1386 } 1387 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1388 } 1389 } 1390 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1391 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1392 } 1393 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1394 } 1395 } 1396 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1397 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1398 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1399 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1400 for (field = 0; field < Nf; ++field) { 1401 PetscObject obj; 1402 PetscClassId id; 1403 PetscDualSpace Q = NULL; 1404 PetscQuadrature f; 1405 const PetscReal *qpoints, *qweights; 1406 PetscInt Nc, Np, fpdim, i, d; 1407 1408 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1409 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1410 if (id == PETSCFE_CLASSID) { 1411 PetscFE fe = (PetscFE) obj; 1412 1413 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1414 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1415 } else if (id == PETSCFV_CLASSID) { 1416 PetscFV fv = (PetscFV) obj; 1417 1418 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1419 Nc = 1; 1420 } 1421 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1422 /* For each fine grid cell */ 1423 for (cell = cStart; cell < cEnd; ++cell) { 1424 PetscInt *findices, *cindices; 1425 PetscInt numFIndices, numCIndices; 1426 1427 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1428 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1429 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1430 for (i = 0; i < fpdim; ++i) { 1431 Vec pointVec; 1432 PetscScalar *pV; 1433 PetscSF coarseCellSF = NULL; 1434 const PetscSFNode *coarseCells; 1435 PetscInt numCoarseCells, cpdim, q, c, j; 1436 1437 /* Get points from the dual basis functional quadrature */ 1438 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1439 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, &qweights);CHKERRQ(ierr); 1440 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1441 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1442 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1443 for (q = 0; q < Np; ++q) { 1444 /* Transform point to real space */ 1445 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1446 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1447 } 1448 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1449 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1450 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1451 /* Update preallocation info */ 1452 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1453 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1454 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1455 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1456 PetscReal pVReal[3]; 1457 1458 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1459 /* Transform points from real space to coarse reference space */ 1460 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 1461 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 1462 CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x); 1463 1464 if (id == PETSCFE_CLASSID) { 1465 PetscFE fe = (PetscFE) obj; 1466 PetscReal *B; 1467 1468 /* Evaluate coarse basis on contained point */ 1469 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1470 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 1471 ierr = PetscMemzero(elemMat, cpdim*Nc*Nc * sizeof(PetscScalar));CHKERRQ(ierr); 1472 /* Get elemMat entries by multiplying by weight */ 1473 for (j = 0; j < cpdim; ++j) { 1474 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = B[j*Nc + c]*qweights[ccell]; 1475 } 1476 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1477 } else { 1478 cpdim = 1; 1479 for (j = 0; j < cpdim; ++j) { 1480 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = 1.0*qweights[ccell]; 1481 } 1482 } 1483 /* Update interpolator */ 1484 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, Nc, numCIndices, elemMat);CHKERRQ(ierr);} 1485 if (numCIndices != cpdim*Nc) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D*%D", numCIndices, cpdim, Nc); 1486 ierr = MatSetValues(In, Nc, &findices[i*Nc], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1487 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1488 } 1489 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1490 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1491 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1492 } 1493 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1494 } 1495 } 1496 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 1497 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 1498 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1499 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1500 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1501 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1502 PetscFunctionReturn(0); 1503 } 1504 1505 /*@ 1506 DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns 1507 1508 Input Parameters: 1509 + dmc - The coarse mesh 1510 - dmf - The fine mesh 1511 - user - The user context 1512 1513 Output Parameter: 1514 . sc - The mapping 1515 1516 Level: developer 1517 1518 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1519 @*/ 1520 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 1521 { 1522 PetscDS prob; 1523 PetscFE *feRef; 1524 PetscFV *fvRef; 1525 Vec fv, cv; 1526 IS fis, cis; 1527 PetscSection fsection, fglobalSection, csection, cglobalSection; 1528 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 1529 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m; 1530 PetscErrorCode ierr; 1531 1532 PetscFunctionBegin; 1533 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1534 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1535 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1536 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1537 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1538 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1539 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1540 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1541 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1542 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1543 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1544 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1545 for (f = 0; f < Nf; ++f) { 1546 PetscObject obj; 1547 PetscClassId id; 1548 PetscInt fNb = 0, Nc = 0; 1549 1550 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1551 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1552 if (id == PETSCFE_CLASSID) { 1553 PetscFE fe = (PetscFE) obj; 1554 1555 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1556 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 1557 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1558 } else if (id == PETSCFV_CLASSID) { 1559 PetscFV fv = (PetscFV) obj; 1560 PetscDualSpace Q; 1561 1562 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1563 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1564 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 1565 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1566 } 1567 fTotDim += fNb*Nc; 1568 } 1569 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1570 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 1571 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 1572 PetscFE feC; 1573 PetscFV fvC; 1574 PetscDualSpace QF, QC; 1575 PetscInt NcF, NcC, fpdim, cpdim; 1576 1577 if (feRef[field]) { 1578 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 1579 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 1580 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 1581 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 1582 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1583 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 1584 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1585 } else { 1586 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 1587 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 1588 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 1589 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 1590 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1591 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 1592 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1593 } 1594 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); 1595 for (c = 0; c < cpdim; ++c) { 1596 PetscQuadrature cfunc; 1597 const PetscReal *cqpoints; 1598 PetscInt NpC; 1599 PetscBool found = PETSC_FALSE; 1600 1601 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 1602 ierr = PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 1603 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 1604 for (f = 0; f < fpdim; ++f) { 1605 PetscQuadrature ffunc; 1606 const PetscReal *fqpoints; 1607 PetscReal sum = 0.0; 1608 PetscInt NpF, comp; 1609 1610 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 1611 ierr = PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 1612 if (NpC != NpF) continue; 1613 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 1614 if (sum > 1.0e-9) continue; 1615 for (comp = 0; comp < NcC; ++comp) { 1616 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 1617 } 1618 found = PETSC_TRUE; 1619 break; 1620 } 1621 if (!found) { 1622 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 1623 if (fvRef[field]) { 1624 PetscInt comp; 1625 for (comp = 0; comp < NcC; ++comp) { 1626 cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp; 1627 } 1628 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 1629 } 1630 } 1631 offsetC += cpdim*NcC; 1632 offsetF += fpdim*NcF; 1633 } 1634 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 1635 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1636 1637 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 1638 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 1639 ierr = VecGetOwnershipRange(cv, &startC, &endC);CHKERRQ(ierr); 1640 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 1641 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 1642 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 1643 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 1644 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 1645 for (c = cStart; c < cEnd; ++c) { 1646 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 1647 for (d = 0; d < cTotDim; ++d) { 1648 if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue; 1649 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]]); 1650 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 1651 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 1652 } 1653 } 1654 ierr = PetscFree(cmap);CHKERRQ(ierr); 1655 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 1656 1657 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 1658 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 1659 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 1660 ierr = ISDestroy(&cis);CHKERRQ(ierr); 1661 ierr = ISDestroy(&fis);CHKERRQ(ierr); 1662 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 1663 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 1664 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1665 PetscFunctionReturn(0); 1666 } 1667