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