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, *v0, *J, *invJ, detJ; 446 PetscReal localDiff = 0.0; 447 const PetscReal *quadPoints, *quadWeights; 448 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 449 PetscErrorCode ierr; 450 451 PetscFunctionBegin; 452 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 453 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 454 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 455 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 456 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 457 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 458 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 459 for (field = 0; field < numFields; ++field) { 460 PetscObject obj; 461 PetscClassId id; 462 PetscInt Nc; 463 464 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 465 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 466 if (id == PETSCFE_CLASSID) { 467 PetscFE fe = (PetscFE) obj; 468 469 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 470 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 471 } else if (id == PETSCFV_CLASSID) { 472 PetscFV fv = (PetscFV) obj; 473 474 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 475 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 476 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 477 numComponents += Nc; 478 } 479 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 480 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 481 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 482 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 483 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 484 for (c = cStart; c < cEnd; ++c) { 485 PetscScalar *x = NULL; 486 PetscReal elemDiff = 0.0; 487 488 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 489 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 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 < numQuadPoints; ++q) { 509 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 510 ierr = (*funcs[field])(dim, time, coords, 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 = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);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);CHKERRQ(ierr);} 523 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 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 = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);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, *v0, *J, *invJ, detJ; 547 PetscReal localDiff = 0.0; 548 PetscInt dim, 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 = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 554 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 555 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 556 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 557 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 558 for (field = 0; field < numFields; ++field) { 559 PetscFE fe; 560 PetscInt Nc; 561 562 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 563 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 564 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 565 numComponents += Nc; 566 } 567 /* ierr = DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 568 ierr = PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr); 569 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 570 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 571 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 572 for (c = cStart; c < cEnd; ++c) { 573 PetscScalar *x = NULL; 574 PetscReal elemDiff = 0.0; 575 576 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 577 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 578 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 579 580 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 581 PetscFE fe; 582 void * const ctx = ctxs ? ctxs[field] : NULL; 583 const PetscReal *quadPoints, *quadWeights; 584 PetscReal *basisDer; 585 PetscInt numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g; 586 587 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 588 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 589 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 590 ierr = PetscFEGetNumComponents(fe, &Ncomp);CHKERRQ(ierr); 591 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 592 if (debug) { 593 char title[1024]; 594 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 595 ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr); 596 } 597 for (q = 0; q < numQuadPoints; ++q) { 598 for (d = 0; d < dim; d++) { 599 coords[d] = v0[d]; 600 for (e = 0; e < dim; e++) { 601 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 602 } 603 } 604 ierr = (*funcs[field])(dim, time, coords, n, numFields, funcVal, ctx); 605 if (ierr) { 606 PetscErrorCode ierr2; 607 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 608 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 609 ierr2 = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr2); 610 CHKERRQ(ierr); 611 } 612 for (fc = 0; fc < Ncomp; ++fc) { 613 PetscScalar interpolant = 0.0; 614 615 for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0; 616 for (f = 0; f < Nb; ++f) { 617 const PetscInt fidx = f*Ncomp+fc; 618 619 for (d = 0; d < dim; ++d) { 620 realSpaceDer[d] = 0.0; 621 for (g = 0; g < dim; ++g) { 622 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g]; 623 } 624 interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d]; 625 } 626 } 627 for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d]; 628 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);CHKERRQ(ierr);} 629 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 630 } 631 } 632 comp += Ncomp; 633 fieldOffset += Nb*Ncomp; 634 } 635 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 636 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 637 localDiff += elemDiff; 638 } 639 ierr = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr); 640 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 641 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 642 *diff = PetscSqrtReal(*diff); 643 PetscFunctionReturn(0); 644 } 645 646 PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 647 { 648 const PetscInt debug = 0; 649 PetscSection section; 650 PetscQuadrature quad; 651 Vec localX; 652 PetscScalar *funcVal, *interpolant; 653 PetscReal *coords, *v0, *J, *invJ, detJ; 654 PetscReal *localDiff; 655 const PetscReal *quadPoints, *quadWeights; 656 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 657 PetscErrorCode ierr; 658 659 PetscFunctionBegin; 660 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 661 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 662 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 663 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 664 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 665 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 666 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 667 for (field = 0; field < numFields; ++field) { 668 PetscObject obj; 669 PetscClassId id; 670 PetscInt Nc; 671 672 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 673 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 674 if (id == PETSCFE_CLASSID) { 675 PetscFE fe = (PetscFE) obj; 676 677 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 678 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 679 } else if (id == PETSCFV_CLASSID) { 680 PetscFV fv = (PetscFV) obj; 681 682 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 683 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 684 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 685 numComponents += Nc; 686 } 687 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 688 ierr = PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 689 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 690 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 691 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 692 for (c = cStart; c < cEnd; ++c) { 693 PetscScalar *x = NULL; 694 695 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 696 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 697 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 698 699 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 700 PetscObject obj; 701 PetscClassId id; 702 void * const ctx = ctxs ? ctxs[field] : NULL; 703 PetscInt Nb, Nc, q, fc; 704 705 PetscReal elemDiff = 0.0; 706 707 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 708 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 709 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 710 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 711 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 712 if (debug) { 713 char title[1024]; 714 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 715 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 716 } 717 for (q = 0; q < numQuadPoints; ++q) { 718 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 719 ierr = (*funcs[field])(dim, time, coords, numFields, funcVal, ctx); 720 if (ierr) { 721 PetscErrorCode ierr2; 722 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 723 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 724 ierr2 = PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr2); 725 CHKERRQ(ierr); 726 } 727 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 728 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 729 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 730 for (fc = 0; fc < Nc; ++fc) { 731 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d field %d point %g %g %g diff %g\n", c, field, dim > 0 ? coords[0] : 0., dim > 1 ? coords[1] : 0., dim > 2 ? coords[2] : 0., PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ);CHKERRQ(ierr);} 732 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 733 } 734 } 735 fieldOffset += Nb*Nc; 736 localDiff[field] += elemDiff; 737 } 738 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 739 } 740 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 741 ierr = MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 742 for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]); 743 ierr = PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 744 PetscFunctionReturn(0); 745 } 746 747 /*@C 748 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. 749 750 Input Parameters: 751 + dm - The DM 752 . time - The time 753 . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation 754 . ctxs - Optional array of contexts to pass to each function, or NULL. 755 - X - The coefficient vector u_h 756 757 Output Parameter: 758 . D - A Vec which holds the difference ||u - u_h||_2 for each cell 759 760 Level: developer 761 762 .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff() 763 @*/ 764 PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D) 765 { 766 PetscSection section; 767 PetscQuadrature quad; 768 Vec localX; 769 PetscScalar *funcVal, *interpolant; 770 PetscReal *coords, *v0, *J, *invJ, detJ; 771 const PetscReal *quadPoints, *quadWeights; 772 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 773 PetscErrorCode ierr; 774 775 PetscFunctionBegin; 776 ierr = VecSet(D, 0.0);CHKERRQ(ierr); 777 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 778 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 779 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 780 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 781 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 782 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 783 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 784 for (field = 0; field < numFields; ++field) { 785 PetscObject obj; 786 PetscClassId id; 787 PetscInt Nc; 788 789 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 790 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 791 if (id == PETSCFE_CLASSID) { 792 PetscFE fe = (PetscFE) obj; 793 794 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 795 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 796 } else if (id == PETSCFV_CLASSID) { 797 PetscFV fv = (PetscFV) obj; 798 799 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 800 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 801 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 802 numComponents += Nc; 803 } 804 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 805 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 806 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 807 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 808 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 809 for (c = cStart; c < cEnd; ++c) { 810 PetscScalar *x = NULL; 811 PetscScalar elemDiff = 0.0; 812 813 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 814 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 815 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 816 817 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 818 PetscObject obj; 819 PetscClassId id; 820 void * const ctx = ctxs ? ctxs[field] : NULL; 821 PetscInt Nb, Nc, q, fc; 822 823 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 824 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 825 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 826 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 827 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 828 if (funcs[field]) { 829 for (q = 0; q < numQuadPoints; ++q) { 830 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 831 ierr = (*funcs[field])(dim, time, coords, Nc, funcVal, ctx); 832 if (ierr) { 833 PetscErrorCode ierr2; 834 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 835 ierr2 = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr2); 836 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 837 CHKERRQ(ierr); 838 } 839 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 840 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 841 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 842 for (fc = 0; fc < Nc; ++fc) { 843 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 844 } 845 } 846 } 847 fieldOffset += Nb*Nc; 848 } 849 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 850 ierr = VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);CHKERRQ(ierr); 851 } 852 ierr = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 853 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 854 ierr = VecSqrtAbs(D);CHKERRQ(ierr); 855 PetscFunctionReturn(0); 856 } 857 858 /*@ 859 DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user 860 861 Input Parameters: 862 + dm - The mesh 863 . X - Local input vector 864 - user - The user context 865 866 Output Parameter: 867 . integral - Local integral for each field 868 869 Level: developer 870 871 .seealso: DMPlexComputeResidualFEM() 872 @*/ 873 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user) 874 { 875 DM_Plex *mesh = (DM_Plex *) dm->data; 876 DM dmAux, dmGrad; 877 Vec localX, A, cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL; 878 PetscDS prob, probAux = NULL; 879 PetscSection section, sectionAux; 880 PetscFV fvm = NULL; 881 PetscFECellGeom *cgeomFEM; 882 PetscFVFaceGeom *fgeomFVM; 883 PetscFVCellGeom *cgeomFVM; 884 PetscScalar *u, *a = NULL; 885 const PetscScalar *lgrad; 886 PetscReal *lintegral; 887 PetscInt *uOff, *uOff_x, *aOff = NULL; 888 PetscInt dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, cEndInterior, c; 889 PetscInt totDim, totDimAux; 890 PetscBool useFVM = PETSC_FALSE; 891 PetscErrorCode ierr; 892 893 PetscFunctionBegin; 894 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 895 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 896 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 897 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 898 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 899 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 900 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 901 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 902 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 903 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 904 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 905 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 906 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 907 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 908 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 909 numCells = cEnd - cStart; 910 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 911 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 912 if (dmAux) { 913 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 914 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 915 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 916 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 917 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 918 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, NULL);CHKERRQ(ierr); 919 } 920 for (f = 0; f < Nf; ++f) { 921 PetscObject obj; 922 PetscClassId id; 923 924 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 925 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 926 if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;} 927 } 928 if (useFVM) { 929 Vec grad; 930 PetscInt fStart, fEnd; 931 PetscBool compGrad; 932 933 ierr = PetscFVGetComputeGradients(fvm, &compGrad);CHKERRQ(ierr); 934 ierr = PetscFVSetComputeGradients(fvm, PETSC_TRUE);CHKERRQ(ierr); 935 ierr = DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);CHKERRQ(ierr); 936 ierr = DMPlexComputeGradientFVM(dm, fvm, faceGeometryFVM, cellGeometryFVM, &dmGrad);CHKERRQ(ierr); 937 ierr = PetscFVSetComputeGradients(fvm, compGrad);CHKERRQ(ierr); 938 ierr = VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);CHKERRQ(ierr); 939 ierr = VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 940 /* Reconstruct and limit cell gradients */ 941 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 942 ierr = DMGetGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 943 ierr = DMPlexReconstructGradients_Internal(dm, fvm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, localX, grad);CHKERRQ(ierr); 944 /* Communicate gradient values */ 945 ierr = DMGetLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 946 ierr = DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 947 ierr = DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 948 ierr = DMRestoreGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 949 /* Handle non-essential (e.g. outflow) boundary values */ 950 ierr = DMPlexInsertBoundaryValues(dm, PETSC_FALSE, localX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);CHKERRQ(ierr); 951 ierr = VecGetArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 952 } 953 ierr = PetscMalloc3(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeomFEM);CHKERRQ(ierr); 954 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 955 for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;} 956 for (c = cStart; c < cEnd; ++c) { 957 PetscScalar *x = NULL; 958 PetscInt i; 959 960 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeomFEM[c].v0, cgeomFEM[c].J, cgeomFEM[c].invJ, &cgeomFEM[c].detJ);CHKERRQ(ierr); 961 if (cgeomFEM[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeomFEM[c].detJ, c); 962 ierr = DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 963 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 964 ierr = DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 965 if (dmAux) { 966 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 967 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 968 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 969 } 970 } 971 for (f = 0; f < Nf; ++f) { 972 PetscObject obj; 973 PetscClassId id; 974 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 975 976 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 977 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 978 if (id == PETSCFE_CLASSID) { 979 PetscFE fe = (PetscFE) obj; 980 PetscQuadrature q; 981 PetscInt Nq, Nb; 982 983 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 984 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 985 ierr = PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 986 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 987 blockSize = Nb*Nq; 988 batchSize = numBlocks * blockSize; 989 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 990 numChunks = numCells / (numBatches*batchSize); 991 Ne = numChunks*numBatches*batchSize; 992 Nr = numCells % (numBatches*batchSize); 993 offset = numCells - Nr; 994 ierr = PetscFEIntegrate(fe, prob, f, Ne, cgeomFEM, u, probAux, a, lintegral);CHKERRQ(ierr); 995 ierr = PetscFEIntegrate(fe, prob, f, Nr, &cgeomFEM[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);CHKERRQ(ierr); 996 } else if (id == PETSCFV_CLASSID) { 997 /* PetscFV fv = (PetscFV) obj; */ 998 PetscInt foff; 999 PetscPointFunc obj_func; 1000 PetscScalar lint; 1001 1002 ierr = PetscDSGetObjective(prob, f, &obj_func);CHKERRQ(ierr); 1003 ierr = PetscDSGetFieldOffset(prob, f, &foff);CHKERRQ(ierr); 1004 if (obj_func) { 1005 for (c = 0; c < numCells; ++c) { 1006 PetscScalar *u_x; 1007 1008 ierr = DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);CHKERRQ(ierr); 1009 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); 1010 lintegral[f] += PetscRealPart(lint)*cgeomFVM[c].volume; 1011 } 1012 } 1013 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 1014 } 1015 if (useFVM) { 1016 ierr = VecRestoreArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1017 ierr = VecRestoreArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);CHKERRQ(ierr); 1018 ierr = VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1019 ierr = DMRestoreLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1020 ierr = VecDestroy(&faceGeometryFVM);CHKERRQ(ierr); 1021 ierr = VecDestroy(&cellGeometryFVM);CHKERRQ(ierr); 1022 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr); 1023 } 1024 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1025 if (mesh->printFEM) { 1026 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");CHKERRQ(ierr); 1027 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);CHKERRQ(ierr);} 1028 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 1029 } 1030 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1031 ierr = MPIU_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr); 1032 ierr = PetscFree3(lintegral,u,cgeomFEM);CHKERRQ(ierr); 1033 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1034 PetscFunctionReturn(0); 1035 } 1036 1037 /*@ 1038 DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1039 1040 Input Parameters: 1041 + dmf - The fine mesh 1042 . dmc - The coarse mesh 1043 - user - The user context 1044 1045 Output Parameter: 1046 . In - The interpolation matrix 1047 1048 Level: developer 1049 1050 .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1051 @*/ 1052 PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user) 1053 { 1054 DM_Plex *mesh = (DM_Plex *) dmc->data; 1055 const char *name = "Interpolator"; 1056 PetscDS prob; 1057 PetscFE *feRef; 1058 PetscFV *fvRef; 1059 PetscSection fsection, fglobalSection; 1060 PetscSection csection, cglobalSection; 1061 PetscScalar *elemMat; 1062 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c; 1063 PetscInt cTotDim, rTotDim = 0; 1064 PetscErrorCode ierr; 1065 1066 PetscFunctionBegin; 1067 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1068 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1069 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1070 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1071 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1072 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1073 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1074 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1075 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1076 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1077 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr); 1078 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1079 for (f = 0; f < Nf; ++f) { 1080 PetscObject obj; 1081 PetscClassId id; 1082 PetscInt rNb = 0, Nc = 0; 1083 1084 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1085 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1086 if (id == PETSCFE_CLASSID) { 1087 PetscFE fe = (PetscFE) obj; 1088 1089 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1090 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1091 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1092 } else if (id == PETSCFV_CLASSID) { 1093 PetscFV fv = (PetscFV) obj; 1094 PetscDualSpace Q; 1095 1096 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1097 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1098 ierr = PetscDualSpaceGetDimension(Q, &rNb);CHKERRQ(ierr); 1099 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1100 } 1101 rTotDim += rNb*Nc; 1102 } 1103 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1104 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 1105 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 1106 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1107 PetscDualSpace Qref; 1108 PetscQuadrature f; 1109 const PetscReal *qpoints, *qweights; 1110 PetscReal *points; 1111 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1112 1113 /* Compose points from all dual basis functionals */ 1114 if (feRef[fieldI]) { 1115 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1116 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 1117 } else { 1118 ierr = PetscFVGetDualSpace(fvRef[fieldI], &Qref);CHKERRQ(ierr); 1119 ierr = PetscFVGetNumComponents(fvRef[fieldI], &Nc);CHKERRQ(ierr); 1120 } 1121 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1122 for (i = 0; i < fpdim; ++i) { 1123 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1124 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1125 npoints += Np; 1126 } 1127 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1128 for (i = 0, k = 0; i < fpdim; ++i) { 1129 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1130 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1131 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1132 } 1133 1134 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1135 PetscObject obj; 1136 PetscClassId id; 1137 PetscReal *B; 1138 PetscInt NcJ = 0, cpdim = 0, j; 1139 1140 ierr = PetscDSGetDiscretization(prob, fieldJ, &obj);CHKERRQ(ierr); 1141 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1142 if (id == PETSCFE_CLASSID) { 1143 PetscFE fe = (PetscFE) obj; 1144 1145 /* Evaluate basis at points */ 1146 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 1147 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1148 /* For now, fields only interpolate themselves */ 1149 if (fieldI == fieldJ) { 1150 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); 1151 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1152 for (i = 0, k = 0; i < fpdim; ++i) { 1153 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1154 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1155 for (p = 0; p < Np; ++p, ++k) { 1156 for (j = 0; j < cpdim; ++j) { 1157 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]; 1158 } 1159 } 1160 } 1161 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1162 } 1163 } else if (id == PETSCFV_CLASSID) { 1164 PetscFV fv = (PetscFV) obj; 1165 1166 /* Evaluate constant function at points */ 1167 ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr); 1168 cpdim = 1; 1169 /* For now, fields only interpolate themselves */ 1170 if (fieldI == fieldJ) { 1171 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); 1172 for (i = 0, k = 0; i < fpdim; ++i) { 1173 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1174 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1175 for (p = 0; p < Np; ++p, ++k) { 1176 for (j = 0; j < cpdim; ++j) { 1177 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p]; 1178 } 1179 } 1180 } 1181 } 1182 } 1183 offsetJ += cpdim*NcJ; 1184 } 1185 offsetI += fpdim*Nc; 1186 ierr = PetscFree(points);CHKERRQ(ierr); 1187 } 1188 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1189 /* Preallocate matrix */ 1190 { 1191 Mat preallocator; 1192 PetscScalar *vals; 1193 PetscInt *cellCIndices, *cellFIndices; 1194 PetscInt locRows, locCols, cell; 1195 1196 ierr = MatGetLocalSize(In, &locRows, &locCols);CHKERRQ(ierr); 1197 ierr = MatCreate(PetscObjectComm((PetscObject) In), &preallocator);CHKERRQ(ierr); 1198 ierr = MatSetType(preallocator, MATPREALLOCATOR);CHKERRQ(ierr); 1199 ierr = MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1200 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 1201 ierr = PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1202 for (cell = cStart; cell < cEnd; ++cell) { 1203 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1204 ierr = MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);CHKERRQ(ierr); 1205 } 1206 ierr = PetscFree3(vals,cellCIndices,cellFIndices);CHKERRQ(ierr); 1207 ierr = MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1208 ierr = MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1209 ierr = MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);CHKERRQ(ierr); 1210 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 1211 } 1212 /* Fill matrix */ 1213 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1214 for (c = cStart; c < cEnd; ++c) { 1215 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1216 } 1217 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1218 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1219 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1220 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1221 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1222 if (mesh->printFEM) { 1223 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1224 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1225 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1226 } 1227 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1228 PetscFunctionReturn(0); 1229 } 1230 1231 /*@ 1232 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 1233 1234 Input Parameters: 1235 + dmf - The fine mesh 1236 . dmc - The coarse mesh 1237 - user - The user context 1238 1239 Output Parameter: 1240 . In - The interpolation matrix 1241 1242 Level: developer 1243 1244 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1245 @*/ 1246 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 1247 { 1248 DM_Plex *mesh = (DM_Plex *) dmf->data; 1249 const char *name = "Interpolator"; 1250 PetscDS prob; 1251 PetscSection fsection, csection, globalFSection, globalCSection; 1252 PetscHashJK ht; 1253 PetscLayout rLayout; 1254 PetscInt *dnz, *onz; 1255 PetscInt locRows, rStart, rEnd; 1256 PetscReal *x, *v0, *J, *invJ, detJ; 1257 PetscReal *v0c, *Jc, *invJc, detJc; 1258 PetscScalar *elemMat; 1259 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1260 PetscErrorCode ierr; 1261 1262 PetscFunctionBegin; 1263 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1264 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1265 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1266 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1267 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1268 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1269 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1270 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1271 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1272 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1273 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1274 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1275 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1276 ierr = PetscMalloc1(totDim*totDim, &elemMat);CHKERRQ(ierr); 1277 1278 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1279 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1280 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1281 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1282 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1283 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1284 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1285 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1286 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1287 for (field = 0; field < Nf; ++field) { 1288 PetscObject obj; 1289 PetscClassId id; 1290 PetscDualSpace Q = NULL; 1291 PetscQuadrature f; 1292 const PetscReal *qpoints; 1293 PetscInt Nc, Np, fpdim, i, d; 1294 1295 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1296 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1297 if (id == PETSCFE_CLASSID) { 1298 PetscFE fe = (PetscFE) obj; 1299 1300 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1301 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1302 } else if (id == PETSCFV_CLASSID) { 1303 PetscFV fv = (PetscFV) obj; 1304 1305 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1306 Nc = 1; 1307 } 1308 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1309 /* For each fine grid cell */ 1310 for (cell = cStart; cell < cEnd; ++cell) { 1311 PetscInt *findices, *cindices; 1312 PetscInt numFIndices, numCIndices; 1313 1314 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1315 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1316 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1317 for (i = 0; i < fpdim; ++i) { 1318 Vec pointVec; 1319 PetscScalar *pV; 1320 PetscSF coarseCellSF = NULL; 1321 const PetscSFNode *coarseCells; 1322 PetscInt numCoarseCells, q, r, c; 1323 1324 /* Get points from the dual basis functional quadrature */ 1325 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1326 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1327 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1328 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1329 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1330 for (q = 0; q < Np; ++q) { 1331 /* Transform point to real space */ 1332 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1333 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1334 } 1335 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1336 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1337 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1338 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 1339 /* Update preallocation info */ 1340 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1341 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1342 for (r = 0; r < Nc; ++r) { 1343 PetscHashJKKey key; 1344 PetscHashJKIter missing, iter; 1345 1346 key.j = findices[i*Nc+r]; 1347 if (key.j < 0) continue; 1348 /* Get indices for coarse elements */ 1349 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1350 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1351 for (c = 0; c < numCIndices; ++c) { 1352 key.k = cindices[c]; 1353 if (key.k < 0) continue; 1354 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1355 if (missing) { 1356 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1357 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1358 else ++onz[key.j-rStart]; 1359 } 1360 } 1361 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1362 } 1363 } 1364 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1365 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1366 } 1367 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1368 } 1369 } 1370 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1371 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1372 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1373 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1374 for (field = 0; field < Nf; ++field) { 1375 PetscObject obj; 1376 PetscClassId id; 1377 PetscDualSpace Q = NULL; 1378 PetscQuadrature f; 1379 const PetscReal *qpoints, *qweights; 1380 PetscInt Nc, Np, fpdim, i, d; 1381 1382 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1383 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1384 if (id == PETSCFE_CLASSID) { 1385 PetscFE fe = (PetscFE) obj; 1386 1387 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1388 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1389 } else if (id == PETSCFV_CLASSID) { 1390 PetscFV fv = (PetscFV) obj; 1391 1392 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1393 Nc = 1; 1394 } 1395 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1396 /* For each fine grid cell */ 1397 for (cell = cStart; cell < cEnd; ++cell) { 1398 PetscInt *findices, *cindices; 1399 PetscInt numFIndices, numCIndices; 1400 1401 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1402 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1403 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1404 for (i = 0; i < fpdim; ++i) { 1405 Vec pointVec; 1406 PetscScalar *pV; 1407 PetscSF coarseCellSF = NULL; 1408 const PetscSFNode *coarseCells; 1409 PetscInt numCoarseCells, cpdim, q, c, j; 1410 1411 /* Get points from the dual basis functional quadrature */ 1412 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1413 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, &qweights);CHKERRQ(ierr); 1414 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1415 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1416 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1417 for (q = 0; q < Np; ++q) { 1418 /* Transform point to real space */ 1419 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1420 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1421 } 1422 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1423 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1424 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1425 /* Update preallocation info */ 1426 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1427 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1428 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1429 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1430 PetscReal pVReal[3]; 1431 1432 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1433 /* Transform points from real space to coarse reference space */ 1434 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 1435 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 1436 CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x); 1437 1438 if (id == PETSCFE_CLASSID) { 1439 PetscFE fe = (PetscFE) obj; 1440 PetscReal *B; 1441 1442 /* Evaluate coarse basis on contained point */ 1443 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1444 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 1445 ierr = PetscMemzero(elemMat, cpdim*Nc*Nc * sizeof(PetscScalar));CHKERRQ(ierr); 1446 /* Get elemMat entries by multiplying by weight */ 1447 for (j = 0; j < cpdim; ++j) { 1448 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = B[j*Nc + c]*qweights[ccell]; 1449 } 1450 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1451 } else { 1452 cpdim = 1; 1453 for (j = 0; j < cpdim; ++j) { 1454 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = 1.0*qweights[ccell]; 1455 } 1456 } 1457 /* Update interpolator */ 1458 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, Nc, numCIndices, elemMat);CHKERRQ(ierr);} 1459 if (numCIndices != cpdim*Nc) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D*%D", numCIndices, cpdim, Nc); 1460 ierr = MatSetValues(In, Nc, &findices[i*Nc], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1461 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1462 } 1463 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1464 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1465 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1466 } 1467 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1468 } 1469 } 1470 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 1471 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 1472 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1473 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1474 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1475 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1476 PetscFunctionReturn(0); 1477 } 1478 1479 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 1480 { 1481 PetscDS prob; 1482 PetscFE *feRef; 1483 PetscFV *fvRef; 1484 Vec fv, cv; 1485 IS fis, cis; 1486 PetscSection fsection, fglobalSection, csection, cglobalSection; 1487 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 1488 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m; 1489 PetscErrorCode ierr; 1490 1491 PetscFunctionBegin; 1492 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1493 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1494 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1495 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1496 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1497 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1498 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1499 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1500 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1501 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1502 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1503 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1504 for (f = 0; f < Nf; ++f) { 1505 PetscObject obj; 1506 PetscClassId id; 1507 PetscInt fNb = 0, Nc = 0; 1508 1509 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1510 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1511 if (id == PETSCFE_CLASSID) { 1512 PetscFE fe = (PetscFE) obj; 1513 1514 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1515 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 1516 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1517 } else if (id == PETSCFV_CLASSID) { 1518 PetscFV fv = (PetscFV) obj; 1519 PetscDualSpace Q; 1520 1521 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1522 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1523 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 1524 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1525 } 1526 fTotDim += fNb*Nc; 1527 } 1528 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1529 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 1530 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 1531 PetscFE feC; 1532 PetscFV fvC; 1533 PetscDualSpace QF, QC; 1534 PetscInt NcF, NcC, fpdim, cpdim; 1535 1536 if (feRef[field]) { 1537 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 1538 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 1539 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 1540 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 1541 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1542 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 1543 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1544 } else { 1545 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 1546 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 1547 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 1548 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 1549 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1550 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 1551 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1552 } 1553 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); 1554 for (c = 0; c < cpdim; ++c) { 1555 PetscQuadrature cfunc; 1556 const PetscReal *cqpoints; 1557 PetscInt NpC; 1558 PetscBool found = PETSC_FALSE; 1559 1560 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 1561 ierr = PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 1562 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 1563 for (f = 0; f < fpdim; ++f) { 1564 PetscQuadrature ffunc; 1565 const PetscReal *fqpoints; 1566 PetscReal sum = 0.0; 1567 PetscInt NpF, comp; 1568 1569 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 1570 ierr = PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 1571 if (NpC != NpF) continue; 1572 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 1573 if (sum > 1.0e-9) continue; 1574 for (comp = 0; comp < NcC; ++comp) { 1575 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 1576 } 1577 found = PETSC_TRUE; 1578 break; 1579 } 1580 if (!found) { 1581 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 1582 if (fvRef[field]) { 1583 PetscInt comp; 1584 for (comp = 0; comp < NcC; ++comp) { 1585 cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp; 1586 } 1587 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 1588 } 1589 } 1590 offsetC += cpdim*NcC; 1591 offsetF += fpdim*NcF; 1592 } 1593 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 1594 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1595 1596 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 1597 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 1598 ierr = VecGetOwnershipRange(cv, &startC, &endC);CHKERRQ(ierr); 1599 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 1600 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 1601 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 1602 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 1603 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 1604 for (c = cStart; c < cEnd; ++c) { 1605 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 1606 for (d = 0; d < cTotDim; ++d) { 1607 if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue; 1608 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]]); 1609 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 1610 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 1611 } 1612 } 1613 ierr = PetscFree(cmap);CHKERRQ(ierr); 1614 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 1615 1616 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 1617 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 1618 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 1619 ierr = ISDestroy(&cis);CHKERRQ(ierr); 1620 ierr = ISDestroy(&fis);CHKERRQ(ierr); 1621 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 1622 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 1623 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1624 PetscFunctionReturn(0); 1625 } 1626