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