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