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