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