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