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