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