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