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