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