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