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