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