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__ "DMPlexComputeIntegralFEM" 1168 /*@ 1169 DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user 1170 1171 Input Parameters: 1172 + dm - The mesh 1173 . X - Local input vector 1174 - user - The user context 1175 1176 Output Parameter: 1177 . integral - Local integral for each field 1178 1179 Level: developer 1180 1181 .seealso: DMPlexComputeResidualFEM() 1182 @*/ 1183 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user) 1184 { 1185 DM_Plex *mesh = (DM_Plex *) dm->data; 1186 DM dmAux, dmGrad; 1187 Vec localX, A, cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL; 1188 PetscDS prob, probAux = NULL; 1189 PetscSection section, sectionAux; 1190 PetscFV fvm = NULL; 1191 PetscFECellGeom *cgeomFEM; 1192 PetscFVFaceGeom *fgeomFVM; 1193 PetscFVCellGeom *cgeomFVM; 1194 PetscScalar *u, *a = NULL; 1195 const PetscScalar *lgrad; 1196 PetscReal *lintegral; 1197 PetscInt *uOff, *uOff_x, *aOff = NULL; 1198 PetscInt dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, cEndInterior, c; 1199 PetscInt totDim, totDimAux; 1200 PetscBool useFVM = PETSC_FALSE; 1201 PetscErrorCode ierr; 1202 1203 PetscFunctionBegin; 1204 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1205 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1206 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 1207 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1208 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1209 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1210 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1211 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1212 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1213 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 1214 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 1215 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1216 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1217 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1218 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1219 numCells = cEnd - cStart; 1220 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1221 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1222 if (dmAux) { 1223 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1224 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 1225 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1226 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1227 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 1228 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, NULL);CHKERRQ(ierr); 1229 } 1230 for (f = 0; f < Nf; ++f) { 1231 PetscObject obj; 1232 PetscClassId id; 1233 1234 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1235 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1236 if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;} 1237 } 1238 if (useFVM) { 1239 Vec grad; 1240 PetscInt fStart, fEnd; 1241 PetscBool compGrad; 1242 1243 ierr = PetscFVGetComputeGradients(fvm, &compGrad);CHKERRQ(ierr); 1244 ierr = PetscFVSetComputeGradients(fvm, PETSC_TRUE);CHKERRQ(ierr); 1245 ierr = DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);CHKERRQ(ierr); 1246 ierr = DMPlexComputeGradientFVM(dm, fvm, faceGeometryFVM, cellGeometryFVM, &dmGrad);CHKERRQ(ierr); 1247 ierr = PetscFVSetComputeGradients(fvm, compGrad);CHKERRQ(ierr); 1248 ierr = VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);CHKERRQ(ierr); 1249 ierr = VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1250 /* Reconstruct and limit cell gradients */ 1251 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1252 ierr = DMGetGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1253 ierr = DMPlexReconstructGradients_Internal(dm, fvm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, localX, grad);CHKERRQ(ierr); 1254 /* Communicate gradient values */ 1255 ierr = DMGetLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1256 ierr = DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1257 ierr = DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1258 ierr = DMRestoreGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1259 /* Handle non-essential (e.g. outflow) boundary values */ 1260 ierr = DMPlexInsertBoundaryValues(dm, PETSC_FALSE, localX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);CHKERRQ(ierr); 1261 ierr = VecGetArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1262 } 1263 ierr = PetscMalloc3(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeomFEM);CHKERRQ(ierr); 1264 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1265 for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;} 1266 for (c = cStart; c < cEnd; ++c) { 1267 PetscScalar *x = NULL; 1268 PetscInt i; 1269 1270 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeomFEM[c].v0, cgeomFEM[c].J, cgeomFEM[c].invJ, &cgeomFEM[c].detJ);CHKERRQ(ierr); 1271 if (cgeomFEM[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeomFEM[c].detJ, c); 1272 ierr = DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1273 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 1274 ierr = DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1275 if (dmAux) { 1276 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1277 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 1278 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1279 } 1280 } 1281 for (f = 0; f < Nf; ++f) { 1282 PetscObject obj; 1283 PetscClassId id; 1284 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 1285 1286 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1287 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1288 if (id == PETSCFE_CLASSID) { 1289 PetscFE fe = (PetscFE) obj; 1290 PetscQuadrature q; 1291 PetscInt Nq, Nb; 1292 1293 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1294 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1295 ierr = PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 1296 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1297 blockSize = Nb*Nq; 1298 batchSize = numBlocks * blockSize; 1299 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1300 numChunks = numCells / (numBatches*batchSize); 1301 Ne = numChunks*numBatches*batchSize; 1302 Nr = numCells % (numBatches*batchSize); 1303 offset = numCells - Nr; 1304 ierr = PetscFEIntegrate(fe, prob, f, Ne, cgeomFEM, u, probAux, a, lintegral);CHKERRQ(ierr); 1305 ierr = PetscFEIntegrate(fe, prob, f, Nr, &cgeomFEM[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);CHKERRQ(ierr); 1306 } else if (id == PETSCFV_CLASSID) { 1307 /* PetscFV fv = (PetscFV) obj; */ 1308 PetscInt foff; 1309 PetscPointFunc obj_func; 1310 PetscScalar lint; 1311 1312 ierr = PetscDSGetObjective(prob, f, &obj_func);CHKERRQ(ierr); 1313 ierr = PetscDSGetFieldOffset(prob, f, &foff);CHKERRQ(ierr); 1314 if (obj_func) { 1315 for (c = 0; c < numCells; ++c) { 1316 PetscScalar *u_x; 1317 1318 ierr = DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);CHKERRQ(ierr); 1319 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); 1320 lintegral[f] += PetscRealPart(lint)*cgeomFVM[c].volume; 1321 } 1322 } 1323 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 1324 } 1325 if (useFVM) { 1326 ierr = VecRestoreArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1327 ierr = VecRestoreArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);CHKERRQ(ierr); 1328 ierr = VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1329 ierr = DMRestoreLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1330 ierr = VecDestroy(&faceGeometryFVM);CHKERRQ(ierr); 1331 ierr = VecDestroy(&cellGeometryFVM);CHKERRQ(ierr); 1332 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr); 1333 } 1334 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1335 if (mesh->printFEM) { 1336 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");CHKERRQ(ierr); 1337 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);CHKERRQ(ierr);} 1338 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 1339 } 1340 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1341 ierr = MPIU_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr); 1342 ierr = PetscFree3(lintegral,u,cgeomFEM);CHKERRQ(ierr); 1343 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1344 PetscFunctionReturn(0); 1345 } 1346 1347 #undef __FUNCT__ 1348 #define __FUNCT__ "DMPlexComputeInterpolatorNested" 1349 /*@ 1350 DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1351 1352 Input Parameters: 1353 + dmf - The fine mesh 1354 . dmc - The coarse mesh 1355 - user - The user context 1356 1357 Output Parameter: 1358 . In - The interpolation matrix 1359 1360 Level: developer 1361 1362 .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1363 @*/ 1364 PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user) 1365 { 1366 DM_Plex *mesh = (DM_Plex *) dmc->data; 1367 const char *name = "Interpolator"; 1368 PetscDS prob; 1369 PetscFE *feRef; 1370 PetscFV *fvRef; 1371 PetscSection fsection, fglobalSection; 1372 PetscSection csection, cglobalSection; 1373 PetscScalar *elemMat; 1374 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c; 1375 PetscInt cTotDim, rTotDim = 0; 1376 PetscErrorCode ierr; 1377 1378 PetscFunctionBegin; 1379 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1380 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1381 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1382 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1383 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1384 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1385 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1386 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1387 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1388 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1389 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr); 1390 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1391 for (f = 0; f < Nf; ++f) { 1392 PetscObject obj; 1393 PetscClassId id; 1394 PetscInt rNb = 0, Nc = 0; 1395 1396 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1397 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1398 if (id == PETSCFE_CLASSID) { 1399 PetscFE fe = (PetscFE) obj; 1400 1401 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1402 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1403 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1404 } else if (id == PETSCFV_CLASSID) { 1405 PetscFV fv = (PetscFV) obj; 1406 PetscDualSpace Q; 1407 1408 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1409 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1410 ierr = PetscDualSpaceGetDimension(Q, &rNb);CHKERRQ(ierr); 1411 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1412 } 1413 rTotDim += rNb*Nc; 1414 } 1415 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1416 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 1417 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 1418 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1419 PetscDualSpace Qref; 1420 PetscQuadrature f; 1421 const PetscReal *qpoints, *qweights; 1422 PetscReal *points; 1423 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1424 1425 /* Compose points from all dual basis functionals */ 1426 if (feRef[fieldI]) { 1427 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1428 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 1429 } else { 1430 ierr = PetscFVGetDualSpace(fvRef[fieldI], &Qref);CHKERRQ(ierr); 1431 ierr = PetscFVGetNumComponents(fvRef[fieldI], &Nc);CHKERRQ(ierr); 1432 } 1433 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1434 for (i = 0; i < fpdim; ++i) { 1435 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1436 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1437 npoints += Np; 1438 } 1439 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1440 for (i = 0, k = 0; i < fpdim; ++i) { 1441 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1442 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1443 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1444 } 1445 1446 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1447 PetscObject obj; 1448 PetscClassId id; 1449 PetscReal *B; 1450 PetscInt NcJ = 0, cpdim = 0, j; 1451 1452 ierr = PetscDSGetDiscretization(prob, fieldJ, &obj);CHKERRQ(ierr); 1453 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1454 if (id == PETSCFE_CLASSID) { 1455 PetscFE fe = (PetscFE) obj; 1456 1457 /* Evaluate basis at points */ 1458 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 1459 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1460 /* For now, fields only interpolate themselves */ 1461 if (fieldI == fieldJ) { 1462 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); 1463 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1464 for (i = 0, k = 0; i < fpdim; ++i) { 1465 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1466 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1467 for (p = 0; p < Np; ++p, ++k) { 1468 for (j = 0; j < cpdim; ++j) { 1469 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]; 1470 } 1471 } 1472 } 1473 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1474 } 1475 } else if (id == PETSCFV_CLASSID) { 1476 PetscFV fv = (PetscFV) obj; 1477 1478 /* Evaluate constant function at points */ 1479 ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr); 1480 cpdim = 1; 1481 /* For now, fields only interpolate themselves */ 1482 if (fieldI == fieldJ) { 1483 if (Nc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in fine space field %d does not match coarse field %d", Nc, NcJ); 1484 for (i = 0, k = 0; i < fpdim; ++i) { 1485 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1486 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1487 for (p = 0; p < Np; ++p, ++k) { 1488 for (j = 0; j < cpdim; ++j) { 1489 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p]; 1490 } 1491 } 1492 } 1493 } 1494 } 1495 offsetJ += cpdim*NcJ; 1496 } 1497 offsetI += fpdim*Nc; 1498 ierr = PetscFree(points);CHKERRQ(ierr); 1499 } 1500 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1501 /* Preallocate matrix */ 1502 { 1503 Mat preallocator; 1504 PetscScalar *vals; 1505 PetscInt *cellCIndices, *cellFIndices; 1506 PetscInt locRows, locCols, cell; 1507 1508 ierr = MatGetLocalSize(In, &locRows, &locCols);CHKERRQ(ierr); 1509 ierr = MatCreate(PetscObjectComm((PetscObject) In), &preallocator);CHKERRQ(ierr); 1510 ierr = MatSetType(preallocator, MATPREALLOCATOR);CHKERRQ(ierr); 1511 ierr = MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1512 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 1513 ierr = PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1514 for (cell = cStart; cell < cEnd; ++cell) { 1515 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1516 ierr = MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);CHKERRQ(ierr); 1517 } 1518 ierr = PetscFree3(vals,cellCIndices,cellFIndices);CHKERRQ(ierr); 1519 ierr = MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1520 ierr = MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1521 ierr = MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);CHKERRQ(ierr); 1522 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 1523 } 1524 /* Fill matrix */ 1525 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1526 for (c = cStart; c < cEnd; ++c) { 1527 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1528 } 1529 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1530 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1531 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1532 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1533 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1534 if (mesh->printFEM) { 1535 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1536 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1537 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1538 } 1539 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1540 PetscFunctionReturn(0); 1541 } 1542 1543 #undef __FUNCT__ 1544 #define __FUNCT__ "DMPlexComputeInterpolatorGeneral" 1545 /*@ 1546 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 1547 1548 Input Parameters: 1549 + dmf - The fine mesh 1550 . dmc - The coarse mesh 1551 - user - The user context 1552 1553 Output Parameter: 1554 . In - The interpolation matrix 1555 1556 Level: developer 1557 1558 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1559 @*/ 1560 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 1561 { 1562 DM_Plex *mesh = (DM_Plex *) dmf->data; 1563 const char *name = "Interpolator"; 1564 PetscDS prob; 1565 PetscSection fsection, csection, globalFSection, globalCSection; 1566 PetscHashJK ht; 1567 PetscLayout rLayout; 1568 PetscInt *dnz, *onz; 1569 PetscInt locRows, rStart, rEnd; 1570 PetscReal *x, *v0, *J, *invJ, detJ; 1571 PetscReal *v0c, *Jc, *invJc, detJc; 1572 PetscScalar *elemMat; 1573 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1574 PetscErrorCode ierr; 1575 1576 PetscFunctionBegin; 1577 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1578 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1579 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1580 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1581 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1582 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1583 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1584 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1585 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1586 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1587 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1588 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1589 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1590 ierr = PetscMalloc1(totDim*totDim, &elemMat);CHKERRQ(ierr); 1591 1592 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1593 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1594 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1595 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1596 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1597 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1598 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1599 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1600 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1601 for (field = 0; field < Nf; ++field) { 1602 PetscObject obj; 1603 PetscClassId id; 1604 PetscDualSpace Q = NULL; 1605 PetscQuadrature f; 1606 const PetscReal *qpoints; 1607 PetscInt Nc, Np, fpdim, i, d; 1608 1609 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1610 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1611 if (id == PETSCFE_CLASSID) { 1612 PetscFE fe = (PetscFE) obj; 1613 1614 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1615 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1616 } else if (id == PETSCFV_CLASSID) { 1617 PetscFV fv = (PetscFV) obj; 1618 1619 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1620 Nc = 1; 1621 } 1622 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1623 /* For each fine grid cell */ 1624 for (cell = cStart; cell < cEnd; ++cell) { 1625 PetscInt *findices, *cindices; 1626 PetscInt numFIndices, numCIndices; 1627 1628 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1629 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1630 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1631 for (i = 0; i < fpdim; ++i) { 1632 Vec pointVec; 1633 PetscScalar *pV; 1634 PetscSF coarseCellSF = NULL; 1635 const PetscSFNode *coarseCells; 1636 PetscInt numCoarseCells, q, r, c; 1637 1638 /* Get points from the dual basis functional quadrature */ 1639 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1640 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1641 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1642 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1643 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1644 for (q = 0; q < Np; ++q) { 1645 /* Transform point to real space */ 1646 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1647 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1648 } 1649 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1650 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1651 ierr = DMLocatePoints(dmc, pointVec, &coarseCellSF);CHKERRQ(ierr); 1652 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 1653 /* Update preallocation info */ 1654 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1655 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1656 for (r = 0; r < Nc; ++r) { 1657 PetscHashJKKey key; 1658 PetscHashJKIter missing, iter; 1659 1660 key.j = findices[i*Nc+r]; 1661 if (key.j < 0) continue; 1662 /* Get indices for coarse elements */ 1663 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1664 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1665 for (c = 0; c < numCIndices; ++c) { 1666 key.k = cindices[c]; 1667 if (key.k < 0) continue; 1668 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1669 if (missing) { 1670 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1671 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1672 else ++onz[key.j-rStart]; 1673 } 1674 } 1675 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1676 } 1677 } 1678 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1679 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1680 } 1681 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1682 } 1683 } 1684 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1685 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1686 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1687 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1688 for (field = 0; field < Nf; ++field) { 1689 PetscObject obj; 1690 PetscClassId id; 1691 PetscDualSpace Q = NULL; 1692 PetscQuadrature f; 1693 const PetscReal *qpoints, *qweights; 1694 PetscInt Nc, Np, fpdim, i, d; 1695 1696 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1697 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1698 if (id == PETSCFE_CLASSID) { 1699 PetscFE fe = (PetscFE) obj; 1700 1701 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1702 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1703 } else if (id == PETSCFV_CLASSID) { 1704 PetscFV fv = (PetscFV) obj; 1705 1706 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1707 Nc = 1; 1708 } 1709 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1710 /* For each fine grid cell */ 1711 for (cell = cStart; cell < cEnd; ++cell) { 1712 PetscInt *findices, *cindices; 1713 PetscInt numFIndices, numCIndices; 1714 1715 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1716 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1717 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1718 for (i = 0; i < fpdim; ++i) { 1719 Vec pointVec; 1720 PetscScalar *pV; 1721 PetscSF coarseCellSF = NULL; 1722 const PetscSFNode *coarseCells; 1723 PetscInt numCoarseCells, cpdim, q, c, j; 1724 1725 /* Get points from the dual basis functional quadrature */ 1726 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1727 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, &qweights);CHKERRQ(ierr); 1728 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1729 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1730 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1731 for (q = 0; q < Np; ++q) { 1732 /* Transform point to real space */ 1733 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1734 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1735 } 1736 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1737 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1738 ierr = DMLocatePoints(dmc, pointVec, &coarseCellSF);CHKERRQ(ierr); 1739 /* Update preallocation info */ 1740 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1741 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1742 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1743 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1744 PetscReal pVReal[3]; 1745 1746 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1747 /* Transform points from real space to coarse reference space */ 1748 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 1749 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 1750 CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x); 1751 1752 if (id == PETSCFE_CLASSID) { 1753 PetscFE fe = (PetscFE) obj; 1754 PetscReal *B; 1755 1756 /* Evaluate coarse basis on contained point */ 1757 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1758 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 1759 /* Get elemMat entries by multiplying by weight */ 1760 for (j = 0; j < cpdim; ++j) { 1761 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = B[j*Nc + c]*qweights[ccell]; 1762 } 1763 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1764 } else { 1765 cpdim = 1; 1766 for (j = 0; j < cpdim; ++j) { 1767 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = 1.0*qweights[ccell]; 1768 } 1769 } 1770 /* Update interpolator */ 1771 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, Nc, numCIndices, elemMat);CHKERRQ(ierr);} 1772 ierr = MatSetValues(In, Nc, &findices[i*Nc], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1773 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1774 } 1775 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1776 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1777 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1778 } 1779 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1780 } 1781 } 1782 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 1783 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 1784 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1785 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1786 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1787 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1788 PetscFunctionReturn(0); 1789 } 1790 1791 #undef __FUNCT__ 1792 #define __FUNCT__ "DMPlexComputeInjectorFEM" 1793 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 1794 { 1795 PetscDS prob; 1796 PetscFE *feRef; 1797 PetscFV *fvRef; 1798 Vec fv, cv; 1799 IS fis, cis; 1800 PetscSection fsection, fglobalSection, csection, cglobalSection; 1801 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 1802 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m; 1803 PetscErrorCode ierr; 1804 1805 PetscFunctionBegin; 1806 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1807 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1808 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1809 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1810 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1811 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1812 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1813 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1814 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1815 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1816 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1817 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1818 for (f = 0; f < Nf; ++f) { 1819 PetscObject obj; 1820 PetscClassId id; 1821 PetscInt fNb = 0, Nc = 0; 1822 1823 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1824 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1825 if (id == PETSCFE_CLASSID) { 1826 PetscFE fe = (PetscFE) obj; 1827 1828 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1829 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 1830 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1831 } else if (id == PETSCFV_CLASSID) { 1832 PetscFV fv = (PetscFV) obj; 1833 PetscDualSpace Q; 1834 1835 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1836 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1837 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 1838 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1839 } 1840 fTotDim += fNb*Nc; 1841 } 1842 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1843 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 1844 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 1845 PetscFE feC; 1846 PetscFV fvC; 1847 PetscDualSpace QF, QC; 1848 PetscInt NcF, NcC, fpdim, cpdim; 1849 1850 if (feRef[field]) { 1851 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 1852 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 1853 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 1854 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 1855 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1856 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 1857 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1858 } else { 1859 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 1860 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 1861 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 1862 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 1863 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1864 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 1865 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1866 } 1867 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); 1868 for (c = 0; c < cpdim; ++c) { 1869 PetscQuadrature cfunc; 1870 const PetscReal *cqpoints; 1871 PetscInt NpC; 1872 PetscBool found = PETSC_FALSE; 1873 1874 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 1875 ierr = PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 1876 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 1877 for (f = 0; f < fpdim; ++f) { 1878 PetscQuadrature ffunc; 1879 const PetscReal *fqpoints; 1880 PetscReal sum = 0.0; 1881 PetscInt NpF, comp; 1882 1883 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 1884 ierr = PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 1885 if (NpC != NpF) continue; 1886 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 1887 if (sum > 1.0e-9) continue; 1888 for (comp = 0; comp < NcC; ++comp) { 1889 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 1890 } 1891 found = PETSC_TRUE; 1892 break; 1893 } 1894 if (!found) { 1895 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 1896 if (fvRef[field]) { 1897 PetscInt comp; 1898 for (comp = 0; comp < NcC; ++comp) { 1899 cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp; 1900 } 1901 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 1902 } 1903 } 1904 offsetC += cpdim*NcC; 1905 offsetF += fpdim*NcF; 1906 } 1907 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 1908 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1909 1910 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 1911 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 1912 ierr = VecGetOwnershipRange(cv, &startC, &endC);CHKERRQ(ierr); 1913 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 1914 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 1915 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 1916 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 1917 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 1918 for (c = cStart; c < cEnd; ++c) { 1919 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 1920 for (d = 0; d < cTotDim; ++d) { 1921 if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue; 1922 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]]); 1923 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 1924 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 1925 } 1926 } 1927 ierr = PetscFree(cmap);CHKERRQ(ierr); 1928 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 1929 1930 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 1931 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 1932 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 1933 ierr = ISDestroy(&cis);CHKERRQ(ierr); 1934 ierr = ISDestroy(&fis);CHKERRQ(ierr); 1935 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 1936 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 1937 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1938 PetscFunctionReturn(0); 1939 } 1940