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