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