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