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 = DMPlexProjectFunction(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(), 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, 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 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, time, &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, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, 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, time, &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, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, 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, time, 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, PetscReal time, PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void *ctx, Vec locX) 601 { 602 PetscErrorCode (**funcs)(PetscInt, PetscReal, 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, time, 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, time, field, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, 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 . time - The time 757 . funcs - The functions to evaluate for each field component 758 . ctxs - Optional array of contexts to pass to each function, or NULL. 759 - X - The coefficient vector u_h 760 761 Output Parameter: 762 . diff - The diff ||u - u_h||_2 763 764 Level: developer 765 766 .seealso: DMPlexProjectFunction(), DMPlexComputeL2FieldDiff(), DMPlexComputeL2GradientDiff() 767 @*/ 768 PetscErrorCode DMPlexComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 769 { 770 const PetscInt debug = 0; 771 PetscSection section; 772 PetscQuadrature quad; 773 Vec localX; 774 PetscScalar *funcVal, *interpolant; 775 PetscReal *coords, *v0, *J, *invJ, detJ; 776 PetscReal localDiff = 0.0; 777 const PetscReal *quadPoints, *quadWeights; 778 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 779 PetscErrorCode ierr; 780 781 PetscFunctionBegin; 782 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 783 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 784 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 785 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 786 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 787 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 788 for (field = 0; field < numFields; ++field) { 789 PetscObject obj; 790 PetscClassId id; 791 PetscInt Nc; 792 793 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 794 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 795 if (id == PETSCFE_CLASSID) { 796 PetscFE fe = (PetscFE) obj; 797 798 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 799 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 800 } else if (id == PETSCFV_CLASSID) { 801 PetscFV fv = (PetscFV) obj; 802 803 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 804 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 805 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 806 numComponents += Nc; 807 } 808 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 809 ierr = DMPlexProjectFunctionLocal(dm, 0.0, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 810 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 811 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 812 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 813 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 814 for (c = cStart; c < cEnd; ++c) { 815 PetscScalar *x = NULL; 816 PetscReal elemDiff = 0.0; 817 818 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 819 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 820 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 821 822 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 823 PetscObject obj; 824 PetscClassId id; 825 void * const ctx = ctxs ? ctxs[field] : NULL; 826 PetscInt Nb, Nc, q, fc; 827 828 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 829 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 830 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 831 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 832 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 833 if (debug) { 834 char title[1024]; 835 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 836 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 837 } 838 for (q = 0; q < numQuadPoints; ++q) { 839 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 840 ierr = (*funcs[field])(dim, time, coords, Nc, funcVal, ctx);CHKERRQ(ierr); 841 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 842 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 843 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 844 for (fc = 0; fc < Nc; ++fc) { 845 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);} 846 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 847 } 848 } 849 fieldOffset += Nb*Nc; 850 } 851 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 852 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 853 localDiff += elemDiff; 854 } 855 ierr = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 856 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 857 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 858 *diff = PetscSqrtReal(*diff); 859 PetscFunctionReturn(0); 860 } 861 862 #undef __FUNCT__ 863 #define __FUNCT__ "DMPlexComputeL2GradientDiff" 864 /*@C 865 DMPlexComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h. 866 867 Input Parameters: 868 + dm - The DM 869 , time - The time 870 . funcs - The gradient functions to evaluate for each field component 871 . ctxs - Optional array of contexts to pass to each function, or NULL. 872 . X - The coefficient vector u_h 873 - n - The vector to project along 874 875 Output Parameter: 876 . diff - The diff ||(grad u - grad u_h) . n||_2 877 878 Level: developer 879 880 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff() 881 @*/ 882 PetscErrorCode DMPlexComputeL2GradientDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff) 883 { 884 const PetscInt debug = 0; 885 PetscSection section; 886 PetscQuadrature quad; 887 Vec localX; 888 PetscScalar *funcVal, *interpolantVec; 889 PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ; 890 PetscReal localDiff = 0.0; 891 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset, comp; 892 PetscErrorCode ierr; 893 894 PetscFunctionBegin; 895 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 896 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 897 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 898 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 899 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 900 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 901 for (field = 0; field < numFields; ++field) { 902 PetscFE fe; 903 PetscInt Nc; 904 905 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 906 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 907 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 908 numComponents += Nc; 909 } 910 /* ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 911 ierr = PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr); 912 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 913 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 914 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 915 for (c = cStart; c < cEnd; ++c) { 916 PetscScalar *x = NULL; 917 PetscReal elemDiff = 0.0; 918 919 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 920 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 921 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 922 923 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 924 PetscFE fe; 925 void * const ctx = ctxs ? ctxs[field] : NULL; 926 const PetscReal *quadPoints, *quadWeights; 927 PetscReal *basisDer; 928 PetscInt numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g; 929 930 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 931 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 932 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 933 ierr = PetscFEGetNumComponents(fe, &Ncomp);CHKERRQ(ierr); 934 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 935 if (debug) { 936 char title[1024]; 937 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 938 ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr); 939 } 940 for (q = 0; q < numQuadPoints; ++q) { 941 for (d = 0; d < dim; d++) { 942 coords[d] = v0[d]; 943 for (e = 0; e < dim; e++) { 944 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 945 } 946 } 947 ierr = (*funcs[field])(dim, time, coords, n, numFields, funcVal, ctx);CHKERRQ(ierr); 948 for (fc = 0; fc < Ncomp; ++fc) { 949 PetscScalar interpolant = 0.0; 950 951 for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0; 952 for (f = 0; f < Nb; ++f) { 953 const PetscInt fidx = f*Ncomp+fc; 954 955 for (d = 0; d < dim; ++d) { 956 realSpaceDer[d] = 0.0; 957 for (g = 0; g < dim; ++g) { 958 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g]; 959 } 960 interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d]; 961 } 962 } 963 for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d]; 964 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);} 965 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 966 } 967 } 968 comp += Ncomp; 969 fieldOffset += Nb*Ncomp; 970 } 971 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 972 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 973 localDiff += elemDiff; 974 } 975 ierr = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr); 976 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 977 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 978 *diff = PetscSqrtReal(*diff); 979 PetscFunctionReturn(0); 980 } 981 982 #undef __FUNCT__ 983 #define __FUNCT__ "DMPlexComputeL2FieldDiff" 984 /*@C 985 DMPlexComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components. 986 987 Input Parameters: 988 + dm - The DM 989 . time - The time 990 . funcs - The functions to evaluate for each field component 991 . ctxs - Optional array of contexts to pass to each function, or NULL. 992 - X - The coefficient vector u_h 993 994 Output Parameter: 995 . diff - The array of differences, ||u^f - u^f_h||_2 996 997 Level: developer 998 999 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff(), DMPlexComputeL2GradientDiff() 1000 @*/ 1001 PetscErrorCode DMPlexComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[]) 1002 { 1003 const PetscInt debug = 0; 1004 PetscSection section; 1005 PetscQuadrature quad; 1006 Vec localX; 1007 PetscScalar *funcVal, *interpolant; 1008 PetscReal *coords, *v0, *J, *invJ, detJ; 1009 PetscReal *localDiff; 1010 const PetscReal *quadPoints, *quadWeights; 1011 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 1012 PetscErrorCode ierr; 1013 1014 PetscFunctionBegin; 1015 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1016 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1017 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1018 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1019 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1020 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1021 for (field = 0; field < numFields; ++field) { 1022 PetscObject obj; 1023 PetscClassId id; 1024 PetscInt Nc; 1025 1026 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1027 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1028 if (id == PETSCFE_CLASSID) { 1029 PetscFE fe = (PetscFE) obj; 1030 1031 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1032 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1033 } else if (id == PETSCFV_CLASSID) { 1034 PetscFV fv = (PetscFV) obj; 1035 1036 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1037 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1038 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1039 numComponents += Nc; 1040 } 1041 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 1042 ierr = DMPlexProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 1043 ierr = PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1044 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1045 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1046 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1047 for (c = cStart; c < cEnd; ++c) { 1048 PetscScalar *x = NULL; 1049 1050 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1051 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 1052 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1053 1054 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1055 PetscObject obj; 1056 PetscClassId id; 1057 void * const ctx = ctxs ? ctxs[field] : NULL; 1058 PetscInt Nb, Nc, q, fc; 1059 1060 PetscReal elemDiff = 0.0; 1061 1062 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1063 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1064 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1065 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1066 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1067 if (debug) { 1068 char title[1024]; 1069 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 1070 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 1071 } 1072 for (q = 0; q < numQuadPoints; ++q) { 1073 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 1074 ierr = (*funcs[field])(dim, time, coords, numFields, funcVal, ctx);CHKERRQ(ierr); 1075 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1076 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1077 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1078 for (fc = 0; fc < Nc; ++fc) { 1079 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);} 1080 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 1081 } 1082 } 1083 fieldOffset += Nb*Nc; 1084 localDiff[field] += elemDiff; 1085 } 1086 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1087 } 1088 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1089 ierr = MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1090 for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]); 1091 ierr = PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 1092 PetscFunctionReturn(0); 1093 } 1094 1095 #undef __FUNCT__ 1096 #define __FUNCT__ "DMPlexComputeL2DiffVec" 1097 /*@C 1098 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. 1099 1100 Input Parameters: 1101 + dm - The DM 1102 . time - The time 1103 . funcs - The functions to evaluate for each field component 1104 . ctxs - Optional array of contexts to pass to each function, or NULL. 1105 - X - The coefficient vector u_h 1106 1107 Output Parameter: 1108 . D - A Vec which holds the difference ||u - u_h||_2 for each cell 1109 1110 Level: developer 1111 1112 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMPlexComputeL2GradientDiff() 1113 @*/ 1114 PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D) 1115 { 1116 PetscSection section; 1117 PetscQuadrature quad; 1118 Vec localX; 1119 PetscScalar *funcVal, *interpolant; 1120 PetscReal *coords, *v0, *J, *invJ, detJ; 1121 const PetscReal *quadPoints, *quadWeights; 1122 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset; 1123 PetscErrorCode ierr; 1124 1125 PetscFunctionBegin; 1126 ierr = VecSet(D, 0.0);CHKERRQ(ierr); 1127 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1128 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1129 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1130 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1131 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1132 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1133 for (field = 0; field < numFields; ++field) { 1134 PetscObject obj; 1135 PetscClassId id; 1136 PetscInt Nc; 1137 1138 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1139 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1140 if (id == PETSCFE_CLASSID) { 1141 PetscFE fe = (PetscFE) obj; 1142 1143 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1144 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1145 } else if (id == PETSCFV_CLASSID) { 1146 PetscFV fv = (PetscFV) obj; 1147 1148 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1149 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1150 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1151 numComponents += Nc; 1152 } 1153 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 1154 ierr = DMPlexProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 1155 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1156 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1157 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1158 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1159 for (c = cStart; c < cEnd; ++c) { 1160 PetscScalar *x = NULL; 1161 PetscScalar elemDiff = 0.0; 1162 1163 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1164 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 1165 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1166 1167 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1168 PetscObject obj; 1169 PetscClassId id; 1170 void * const ctx = ctxs ? ctxs[field] : NULL; 1171 PetscInt Nb, Nc, q, fc; 1172 1173 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 1174 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1175 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1176 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1177 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1178 for (q = 0; q < numQuadPoints; ++q) { 1179 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 1180 ierr = (*funcs[field])(dim, time, coords, Nc, funcVal, ctx);CHKERRQ(ierr); 1181 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1182 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1183 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1184 for (fc = 0; fc < Nc; ++fc) { 1185 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 1186 } 1187 } 1188 fieldOffset += Nb*Nc; 1189 } 1190 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1191 ierr = VecSetValuesSection(D, section, c, &elemDiff, ADD_VALUES);CHKERRQ(ierr); 1192 } 1193 ierr = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 1194 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1195 ierr = VecSqrtAbs(D);CHKERRQ(ierr); 1196 PetscFunctionReturn(0); 1197 } 1198 1199 #undef __FUNCT__ 1200 #define __FUNCT__ "DMPlexComputeIntegralFEM" 1201 /*@ 1202 DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user 1203 1204 Input Parameters: 1205 + dm - The mesh 1206 . X - Local input vector 1207 - user - The user context 1208 1209 Output Parameter: 1210 . integral - Local integral for each field 1211 1212 Level: developer 1213 1214 .seealso: DMPlexComputeResidualFEM() 1215 @*/ 1216 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user) 1217 { 1218 DM_Plex *mesh = (DM_Plex *) dm->data; 1219 DM dmAux; 1220 Vec localX, A; 1221 PetscDS prob, probAux = NULL; 1222 PetscSection section, sectionAux; 1223 PetscFECellGeom *cgeom; 1224 PetscScalar *u, *a = NULL; 1225 PetscReal *lintegral, *vol; 1226 PetscInt dim, Nf, f, numCells, cStart, cEnd, cEndInterior, c; 1227 PetscInt totDim, totDimAux; 1228 PetscErrorCode ierr; 1229 1230 PetscFunctionBegin; 1231 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1232 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1233 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1234 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1235 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1236 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1237 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1238 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1239 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1240 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1241 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1242 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1243 numCells = cEnd - cStart; 1244 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1245 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1246 if (dmAux) { 1247 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1248 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1249 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1250 } 1251 ierr = DMPlexInsertBoundaryValues(dm, localX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 1252 ierr = PetscMalloc4(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeom,numCells,&vol);CHKERRQ(ierr); 1253 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1254 for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;} 1255 for (c = cStart; c < cEnd; ++c) { 1256 PetscScalar *x = NULL; 1257 PetscInt i; 1258 1259 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeom[c].v0, cgeom[c].J, cgeom[c].invJ, &cgeom[c].detJ);CHKERRQ(ierr); 1260 ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol[c], NULL, NULL);CHKERRQ(ierr); 1261 if (cgeom[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeom[c].detJ, c); 1262 ierr = DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1263 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 1264 ierr = DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1265 if (dmAux) { 1266 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1267 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 1268 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1269 } 1270 } 1271 for (f = 0; f < Nf; ++f) { 1272 PetscObject obj; 1273 PetscClassId id; 1274 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 1275 1276 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1277 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1278 if (id == PETSCFE_CLASSID) { 1279 PetscFE fe = (PetscFE) obj; 1280 PetscQuadrature q; 1281 PetscInt Nq, Nb; 1282 1283 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1284 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1285 ierr = PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 1286 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1287 blockSize = Nb*Nq; 1288 batchSize = numBlocks * blockSize; 1289 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1290 numChunks = numCells / (numBatches*batchSize); 1291 Ne = numChunks*numBatches*batchSize; 1292 Nr = numCells % (numBatches*batchSize); 1293 offset = numCells - Nr; 1294 ierr = PetscFEIntegrate(fe, prob, f, Ne, cgeom, u, probAux, a, lintegral);CHKERRQ(ierr); 1295 ierr = PetscFEIntegrate(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);CHKERRQ(ierr); 1296 } else if (id == PETSCFV_CLASSID) { 1297 /* PetscFV fv = (PetscFV) obj; */ 1298 PetscInt foff; 1299 PetscPointFunc obj_func; 1300 PetscScalar lint; 1301 1302 ierr = PetscDSGetObjective(prob, f, &obj_func);CHKERRQ(ierr); 1303 ierr = PetscDSGetFieldOffset(prob, f, &foff);CHKERRQ(ierr); 1304 if (obj_func) { 1305 for (c = 0; c < numCells; ++c) { 1306 /* TODO: Need full pointwise interpolation and get centroid for x argument */ 1307 obj_func(dim, Nf, 0, NULL, NULL, &u[totDim*c+foff], NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0.0, NULL, &lint); 1308 lintegral[f] = PetscRealPart(lint)*vol[c]; 1309 } 1310 } 1311 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 1312 } 1313 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1314 if (mesh->printFEM) { 1315 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");CHKERRQ(ierr); 1316 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);CHKERRQ(ierr);} 1317 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 1318 } 1319 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1320 ierr = MPIU_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr); 1321 ierr = PetscFree4(lintegral,u,cgeom,vol);CHKERRQ(ierr); 1322 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1323 PetscFunctionReturn(0); 1324 } 1325 1326 #undef __FUNCT__ 1327 #define __FUNCT__ "DMPlexComputeInterpolatorNested" 1328 /*@ 1329 DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1330 1331 Input Parameters: 1332 + dmf - The fine mesh 1333 . dmc - The coarse mesh 1334 - user - The user context 1335 1336 Output Parameter: 1337 . In - The interpolation matrix 1338 1339 Level: developer 1340 1341 .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1342 @*/ 1343 PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user) 1344 { 1345 DM_Plex *mesh = (DM_Plex *) dmc->data; 1346 const char *name = "Interpolator"; 1347 PetscDS prob; 1348 PetscFE *feRef; 1349 PetscFV *fvRef; 1350 PetscSection fsection, fglobalSection; 1351 PetscSection csection, cglobalSection; 1352 PetscScalar *elemMat; 1353 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c; 1354 PetscInt cTotDim, rTotDim = 0; 1355 PetscErrorCode ierr; 1356 1357 PetscFunctionBegin; 1358 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1359 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1360 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1361 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1362 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1363 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1364 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1365 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1366 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1367 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1368 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr); 1369 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1370 for (f = 0; f < Nf; ++f) { 1371 PetscObject obj; 1372 PetscClassId id; 1373 PetscInt rNb = 0, Nc = 0; 1374 1375 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1376 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1377 if (id == PETSCFE_CLASSID) { 1378 PetscFE fe = (PetscFE) obj; 1379 1380 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1381 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1382 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1383 } else if (id == PETSCFV_CLASSID) { 1384 PetscFV fv = (PetscFV) obj; 1385 PetscDualSpace Q; 1386 1387 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1388 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1389 ierr = PetscDualSpaceGetDimension(Q, &rNb);CHKERRQ(ierr); 1390 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1391 } 1392 rTotDim += rNb*Nc; 1393 } 1394 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1395 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 1396 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 1397 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1398 PetscDualSpace Qref; 1399 PetscQuadrature f; 1400 const PetscReal *qpoints, *qweights; 1401 PetscReal *points; 1402 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1403 1404 /* Compose points from all dual basis functionals */ 1405 if (feRef[fieldI]) { 1406 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1407 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 1408 } else { 1409 ierr = PetscFVGetDualSpace(fvRef[fieldI], &Qref);CHKERRQ(ierr); 1410 ierr = PetscFVGetNumComponents(fvRef[fieldI], &Nc);CHKERRQ(ierr); 1411 } 1412 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1413 for (i = 0; i < fpdim; ++i) { 1414 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1415 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1416 npoints += Np; 1417 } 1418 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1419 for (i = 0, k = 0; i < fpdim; ++i) { 1420 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1421 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1422 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1423 } 1424 1425 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1426 PetscObject obj; 1427 PetscClassId id; 1428 PetscReal *B; 1429 PetscInt NcJ = 0, cpdim = 0, j; 1430 1431 ierr = PetscDSGetDiscretization(prob, fieldJ, &obj);CHKERRQ(ierr); 1432 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1433 if (id == PETSCFE_CLASSID) { 1434 PetscFE fe = (PetscFE) obj; 1435 1436 /* Evaluate basis at points */ 1437 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 1438 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1439 /* For now, fields only interpolate themselves */ 1440 if (fieldI == fieldJ) { 1441 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); 1442 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1443 for (i = 0, k = 0; i < fpdim; ++i) { 1444 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1445 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1446 for (p = 0; p < Np; ++p, ++k) { 1447 for (j = 0; j < cpdim; ++j) { 1448 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]; 1449 } 1450 } 1451 } 1452 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1453 } 1454 } else if (id == PETSCFV_CLASSID) { 1455 PetscFV fv = (PetscFV) obj; 1456 1457 /* Evaluate constant function at points */ 1458 ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr); 1459 cpdim = 1; 1460 /* For now, fields only interpolate themselves */ 1461 if (fieldI == fieldJ) { 1462 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); 1463 for (i = 0, k = 0; i < fpdim; ++i) { 1464 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1465 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1466 for (p = 0; p < Np; ++p, ++k) { 1467 for (j = 0; j < cpdim; ++j) { 1468 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p]; 1469 } 1470 } 1471 } 1472 } 1473 } 1474 offsetJ += cpdim*NcJ; 1475 } 1476 offsetI += fpdim*Nc; 1477 ierr = PetscFree(points);CHKERRQ(ierr); 1478 } 1479 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1480 /* Preallocate matrix */ 1481 { 1482 PetscHashJK ht; 1483 PetscLayout rLayout; 1484 PetscInt *dnz, *onz, *cellCIndices, *cellFIndices; 1485 PetscInt locRows, rStart, rEnd, cell, r; 1486 1487 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1488 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1489 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1490 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1491 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1492 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1493 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1494 ierr = PetscCalloc4(locRows,&dnz,locRows,&onz,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1495 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1496 for (cell = cStart; cell < cEnd; ++cell) { 1497 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1498 for (r = 0; r < rTotDim; ++r) { 1499 PetscHashJKKey key; 1500 PetscHashJKIter missing, iter; 1501 1502 key.j = cellFIndices[r]; 1503 if (key.j < 0) continue; 1504 for (c = 0; c < cTotDim; ++c) { 1505 key.k = cellCIndices[c]; 1506 if (key.k < 0) continue; 1507 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1508 if (missing) { 1509 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1510 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1511 else ++onz[key.j-rStart]; 1512 } 1513 } 1514 } 1515 } 1516 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1517 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1518 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1519 ierr = PetscFree4(dnz,onz,cellCIndices,cellFIndices);CHKERRQ(ierr); 1520 } 1521 /* Fill matrix */ 1522 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1523 for (c = cStart; c < cEnd; ++c) { 1524 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1525 } 1526 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1527 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1528 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1529 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1530 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1531 if (mesh->printFEM) { 1532 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1533 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1534 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1535 } 1536 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1537 PetscFunctionReturn(0); 1538 } 1539 1540 #undef __FUNCT__ 1541 #define __FUNCT__ "DMPlexComputeInterpolatorGeneral" 1542 /*@ 1543 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 1544 1545 Input Parameters: 1546 + dmf - The fine mesh 1547 . dmc - The coarse mesh 1548 - user - The user context 1549 1550 Output Parameter: 1551 . In - The interpolation matrix 1552 1553 Level: developer 1554 1555 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1556 @*/ 1557 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 1558 { 1559 PetscDS prob; 1560 PetscSection fsection, csection, globalFSection, globalCSection; 1561 PetscHashJK ht; 1562 PetscLayout rLayout; 1563 PetscInt *dnz, *onz; 1564 PetscInt locRows, rStart, rEnd; 1565 PetscReal *x, *v0, *J, *invJ, detJ; 1566 PetscReal *v0c, *Jc, *invJc, detJc; 1567 PetscScalar *elemMat; 1568 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1569 PetscErrorCode ierr; 1570 1571 PetscFunctionBegin; 1572 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1573 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1574 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1575 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1576 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1577 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1578 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1579 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1580 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1581 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1582 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1583 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1584 ierr = PetscMalloc1(totDim*totDim, &elemMat);CHKERRQ(ierr); 1585 1586 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1587 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1588 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1589 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1590 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1591 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1592 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1593 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1594 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1595 for (field = 0; field < Nf; ++field) { 1596 PetscObject obj; 1597 PetscClassId id; 1598 PetscDualSpace Q = NULL; 1599 PetscQuadrature f; 1600 const PetscReal *qpoints; 1601 PetscInt Nc, Np, fpdim, i, d; 1602 1603 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1604 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1605 if (id == PETSCFE_CLASSID) { 1606 PetscFE fe = (PetscFE) obj; 1607 1608 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1609 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1610 } else if (id == PETSCFV_CLASSID) { 1611 PetscFV fv = (PetscFV) obj; 1612 1613 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1614 Nc = 1; 1615 } 1616 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1617 /* For each fine grid cell */ 1618 for (cell = cStart; cell < cEnd; ++cell) { 1619 PetscInt *findices, *cindices; 1620 PetscInt numFIndices, numCIndices; 1621 1622 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1623 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1624 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1625 for (i = 0; i < fpdim; ++i) { 1626 Vec pointVec; 1627 PetscScalar *pV; 1628 IS coarseCellIS; 1629 const PetscInt *coarseCells; 1630 PetscInt numCoarseCells, q, r, c; 1631 1632 /* Get points from the dual basis functional quadrature */ 1633 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1634 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1635 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1636 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1637 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1638 for (q = 0; q < Np; ++q) { 1639 /* Transform point to real space */ 1640 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1641 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1642 } 1643 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1644 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1645 ierr = DMLocatePoints(dmc, pointVec, &coarseCellIS);CHKERRQ(ierr); 1646 ierr = ISViewFromOptions(coarseCellIS, NULL, "-interp_is_view");CHKERRQ(ierr); 1647 /* Update preallocation info */ 1648 ierr = ISGetLocalSize(coarseCellIS, &numCoarseCells);CHKERRQ(ierr); 1649 ierr = ISGetIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1650 for (r = 0; r < Nc; ++r) { 1651 PetscHashJKKey key; 1652 PetscHashJKIter missing, iter; 1653 1654 key.j = findices[i*Nc+r]; 1655 if (key.j < 0) continue; 1656 /* Get indices for coarse elements */ 1657 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1658 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1659 for (c = 0; c < numCIndices; ++c) { 1660 key.k = cindices[c]; 1661 if (key.k < 0) continue; 1662 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1663 if (missing) { 1664 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1665 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1666 else ++onz[key.j-rStart]; 1667 } 1668 } 1669 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1670 } 1671 } 1672 ierr = ISRestoreIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1673 ierr = ISDestroy(&coarseCellIS);CHKERRQ(ierr); 1674 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1675 } 1676 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1677 } 1678 } 1679 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1680 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1681 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1682 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1683 for (field = 0; field < Nf; ++field) { 1684 PetscObject obj; 1685 PetscClassId id; 1686 PetscDualSpace Q = NULL; 1687 PetscQuadrature f; 1688 const PetscReal *qpoints, *qweights; 1689 PetscInt Nc, Np, fpdim, i, d; 1690 1691 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1692 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1693 if (id == PETSCFE_CLASSID) { 1694 PetscFE fe = (PetscFE) obj; 1695 1696 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1697 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1698 } else if (id == PETSCFV_CLASSID) { 1699 PetscFV fv = (PetscFV) obj; 1700 1701 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1702 Nc = 1; 1703 } 1704 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1705 /* For each fine grid cell */ 1706 for (cell = cStart; cell < cEnd; ++cell) { 1707 PetscInt *findices, *cindices; 1708 PetscInt numFIndices, numCIndices; 1709 1710 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1711 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1712 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1713 for (i = 0; i < fpdim; ++i) { 1714 Vec pointVec; 1715 PetscScalar *pV; 1716 IS coarseCellIS; 1717 const PetscInt *coarseCells; 1718 PetscInt numCoarseCells, cpdim, q, c, j; 1719 1720 /* Get points from the dual basis functional quadrature */ 1721 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1722 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, &qweights);CHKERRQ(ierr); 1723 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1724 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1725 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1726 for (q = 0; q < Np; ++q) { 1727 /* Transform point to real space */ 1728 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1729 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1730 } 1731 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1732 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1733 ierr = DMLocatePoints(dmc, pointVec, &coarseCellIS);CHKERRQ(ierr); 1734 /* Update preallocation info */ 1735 ierr = ISGetLocalSize(coarseCellIS, &numCoarseCells);CHKERRQ(ierr); 1736 ierr = ISGetIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1737 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1738 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1739 PetscReal pVReal[3]; 1740 1741 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1742 /* Transform points from real space to coarse reference space */ 1743 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell], NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 1744 CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x); 1745 for (d = 0; d < dim; ++d) pV[ccell*dim+d] = pVReal[d]; 1746 1747 if (id == PETSCFE_CLASSID) { 1748 PetscFE fe = (PetscFE) obj; 1749 PetscReal *B; 1750 1751 /* Evaluate coarse basis on contained point */ 1752 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1753 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 1754 /* Get elemMat entries by multiplying by weight */ 1755 for (j = 0; j < cpdim; ++j) { 1756 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = B[j*Nc + c]*qweights[ccell]; 1757 } 1758 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1759 } else { 1760 cpdim = 1; 1761 for (j = 0; j < cpdim; ++j) { 1762 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = 1.0*qweights[ccell]; 1763 } 1764 } 1765 /* Update interpolator */ 1766 ierr = MatSetValues(In, Nc, &findices[i*Nc], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1767 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1768 } 1769 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1770 ierr = ISRestoreIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1771 ierr = ISDestroy(&coarseCellIS);CHKERRQ(ierr); 1772 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1773 } 1774 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1775 } 1776 } 1777 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 1778 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 1779 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1780 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1781 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1782 PetscFunctionReturn(0); 1783 } 1784 1785 #undef __FUNCT__ 1786 #define __FUNCT__ "DMPlexComputeInjectorFEM" 1787 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 1788 { 1789 PetscDS prob; 1790 PetscFE *feRef; 1791 PetscFV *fvRef; 1792 Vec fv, cv; 1793 IS fis, cis; 1794 PetscSection fsection, fglobalSection, csection, cglobalSection; 1795 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 1796 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, offsetC, offsetF, m; 1797 PetscErrorCode ierr; 1798 1799 PetscFunctionBegin; 1800 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1801 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1802 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1803 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1804 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1805 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1806 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1807 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1808 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1809 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1810 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1811 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1812 for (f = 0; f < Nf; ++f) { 1813 PetscObject obj; 1814 PetscClassId id; 1815 PetscInt fNb = 0, Nc = 0; 1816 1817 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1818 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1819 if (id == PETSCFE_CLASSID) { 1820 PetscFE fe = (PetscFE) obj; 1821 1822 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1823 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 1824 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1825 } else if (id == PETSCFV_CLASSID) { 1826 PetscFV fv = (PetscFV) obj; 1827 PetscDualSpace Q; 1828 1829 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1830 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1831 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 1832 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1833 } 1834 fTotDim += fNb*Nc; 1835 } 1836 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1837 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 1838 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 1839 PetscFE feC; 1840 PetscFV fvC; 1841 PetscDualSpace QF, QC; 1842 PetscInt NcF, NcC, fpdim, cpdim; 1843 1844 if (feRef[field]) { 1845 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 1846 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 1847 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 1848 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 1849 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1850 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 1851 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1852 } else { 1853 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 1854 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 1855 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 1856 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 1857 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1858 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 1859 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1860 } 1861 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); 1862 for (c = 0; c < cpdim; ++c) { 1863 PetscQuadrature cfunc; 1864 const PetscReal *cqpoints; 1865 PetscInt NpC; 1866 PetscBool found = PETSC_FALSE; 1867 1868 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 1869 ierr = PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 1870 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 1871 for (f = 0; f < fpdim; ++f) { 1872 PetscQuadrature ffunc; 1873 const PetscReal *fqpoints; 1874 PetscReal sum = 0.0; 1875 PetscInt NpF, comp; 1876 1877 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 1878 ierr = PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 1879 if (NpC != NpF) continue; 1880 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 1881 if (sum > 1.0e-9) continue; 1882 for (comp = 0; comp < NcC; ++comp) { 1883 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 1884 } 1885 found = PETSC_TRUE; 1886 break; 1887 } 1888 if (!found) { 1889 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 1890 if (fvRef[field]) { 1891 PetscInt comp; 1892 for (comp = 0; comp < NcC; ++comp) { 1893 cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp; 1894 } 1895 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 1896 } 1897 } 1898 offsetC += cpdim*NcC; 1899 offsetF += fpdim*NcF; 1900 } 1901 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 1902 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1903 1904 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 1905 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 1906 ierr = VecGetOwnershipRange(cv, &startC, NULL);CHKERRQ(ierr); 1907 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 1908 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 1909 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 1910 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 1911 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 1912 for (c = cStart; c < cEnd; ++c) { 1913 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 1914 for (d = 0; d < cTotDim; ++d) { 1915 if (cellCIndices[d] < 0) continue; 1916 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]]); 1917 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 1918 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 1919 } 1920 } 1921 ierr = PetscFree(cmap);CHKERRQ(ierr); 1922 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 1923 1924 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 1925 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 1926 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 1927 ierr = ISDestroy(&cis);CHKERRQ(ierr); 1928 ierr = ISDestroy(&fis);CHKERRQ(ierr); 1929 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 1930 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 1931 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1932 PetscFunctionReturn(0); 1933 } 1934