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