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