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