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, numFields, 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__ "DMPlexComputeIntegralFEM" 1088 /*@ 1089 DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user 1090 1091 Input Parameters: 1092 + dm - The mesh 1093 . X - Local input vector 1094 - user - The user context 1095 1096 Output Parameter: 1097 . integral - Local integral for each field 1098 1099 Level: developer 1100 1101 .seealso: DMPlexComputeResidualFEM() 1102 @*/ 1103 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user) 1104 { 1105 DM_Plex *mesh = (DM_Plex *) dm->data; 1106 DM dmAux; 1107 Vec localX, A; 1108 PetscDS prob, probAux = NULL; 1109 PetscSection section, sectionAux; 1110 PetscFECellGeom *cgeom; 1111 PetscScalar *u, *a = NULL; 1112 PetscReal *lintegral, *vol; 1113 PetscInt dim, Nf, f, numCells, cStart, cEnd, cEndInterior, c; 1114 PetscInt totDim, totDimAux; 1115 PetscErrorCode ierr; 1116 1117 PetscFunctionBegin; 1118 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1119 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1120 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1121 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1122 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1123 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1124 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1125 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1126 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1127 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1128 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1129 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1130 numCells = cEnd - cStart; 1131 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1132 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1133 if (dmAux) { 1134 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1135 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1136 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1137 } 1138 ierr = DMPlexInsertBoundaryValues(dm, localX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 1139 ierr = PetscMalloc4(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeom,numCells,&vol);CHKERRQ(ierr); 1140 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1141 for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;} 1142 for (c = cStart; c < cEnd; ++c) { 1143 PetscScalar *x = NULL; 1144 PetscInt i; 1145 1146 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeom[c].v0, cgeom[c].J, cgeom[c].invJ, &cgeom[c].detJ);CHKERRQ(ierr); 1147 ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol[c], NULL, NULL);CHKERRQ(ierr); 1148 if (cgeom[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeom[c].detJ, c); 1149 ierr = DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1150 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 1151 ierr = DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1152 if (dmAux) { 1153 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1154 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 1155 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1156 } 1157 } 1158 for (f = 0; f < Nf; ++f) { 1159 PetscObject obj; 1160 PetscClassId id; 1161 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 1162 1163 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1164 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1165 if (id == PETSCFE_CLASSID) { 1166 PetscFE fe = (PetscFE) obj; 1167 PetscQuadrature q; 1168 PetscInt Nq, Nb; 1169 1170 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1171 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1172 ierr = PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 1173 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1174 blockSize = Nb*Nq; 1175 batchSize = numBlocks * blockSize; 1176 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1177 numChunks = numCells / (numBatches*batchSize); 1178 Ne = numChunks*numBatches*batchSize; 1179 Nr = numCells % (numBatches*batchSize); 1180 offset = numCells - Nr; 1181 ierr = PetscFEIntegrate(fe, prob, f, Ne, cgeom, u, probAux, a, lintegral);CHKERRQ(ierr); 1182 ierr = PetscFEIntegrate(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);CHKERRQ(ierr); 1183 } else if (id == PETSCFV_CLASSID) { 1184 /* PetscFV fv = (PetscFV) obj; */ 1185 PetscInt foff; 1186 PetscPointFunc obj_func; 1187 PetscScalar lint; 1188 1189 ierr = PetscDSGetObjective(prob, f, &obj_func);CHKERRQ(ierr); 1190 ierr = PetscDSGetFieldOffset(prob, f, &foff);CHKERRQ(ierr); 1191 if (obj_func) { 1192 for (c = 0; c < numCells; ++c) { 1193 /* TODO: Need full pointwise interpolation and get centroid for x argument */ 1194 obj_func(dim, Nf, 0, NULL, NULL, &u[totDim*c+foff], NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0.0, NULL, &lint); 1195 lintegral[f] = PetscRealPart(lint)*vol[c]; 1196 } 1197 } 1198 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 1199 } 1200 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1201 if (mesh->printFEM) { 1202 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");CHKERRQ(ierr); 1203 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);CHKERRQ(ierr);} 1204 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 1205 } 1206 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1207 ierr = MPI_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr); 1208 ierr = PetscFree4(lintegral,u,cgeom,vol);CHKERRQ(ierr); 1209 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1210 PetscFunctionReturn(0); 1211 } 1212 1213 #undef __FUNCT__ 1214 #define __FUNCT__ "DMPlexComputeInterpolatorNested" 1215 /*@ 1216 DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1217 1218 Input Parameters: 1219 + dmf - The fine mesh 1220 . dmc - The coarse mesh 1221 - user - The user context 1222 1223 Output Parameter: 1224 . In - The interpolation matrix 1225 1226 Level: developer 1227 1228 .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1229 @*/ 1230 PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user) 1231 { 1232 DM_Plex *mesh = (DM_Plex *) dmc->data; 1233 const char *name = "Interpolator"; 1234 PetscDS prob; 1235 PetscFE *feRef; 1236 PetscFV *fvRef; 1237 PetscSection fsection, fglobalSection; 1238 PetscSection csection, cglobalSection; 1239 PetscScalar *elemMat; 1240 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c; 1241 PetscInt cTotDim, rTotDim = 0; 1242 PetscErrorCode ierr; 1243 1244 PetscFunctionBegin; 1245 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1246 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1247 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1248 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1249 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1250 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1251 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1252 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1253 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1254 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1255 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr); 1256 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1257 for (f = 0; f < Nf; ++f) { 1258 PetscObject obj; 1259 PetscClassId id; 1260 PetscInt rNb = 0, Nc = 0; 1261 1262 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1263 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1264 if (id == PETSCFE_CLASSID) { 1265 PetscFE fe = (PetscFE) obj; 1266 1267 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1268 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1269 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1270 } else if (id == PETSCFV_CLASSID) { 1271 PetscFV fv = (PetscFV) obj; 1272 PetscDualSpace Q; 1273 1274 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1275 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1276 ierr = PetscDualSpaceGetDimension(Q, &rNb);CHKERRQ(ierr); 1277 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1278 } 1279 rTotDim += rNb*Nc; 1280 } 1281 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1282 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 1283 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 1284 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1285 PetscDualSpace Qref; 1286 PetscQuadrature f; 1287 const PetscReal *qpoints, *qweights; 1288 PetscReal *points; 1289 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1290 1291 /* Compose points from all dual basis functionals */ 1292 if (feRef[fieldI]) { 1293 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1294 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 1295 } else { 1296 ierr = PetscFVGetDualSpace(fvRef[fieldI], &Qref);CHKERRQ(ierr); 1297 ierr = PetscFVGetNumComponents(fvRef[fieldI], &Nc);CHKERRQ(ierr); 1298 } 1299 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1300 for (i = 0; i < fpdim; ++i) { 1301 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1302 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1303 npoints += Np; 1304 } 1305 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1306 for (i = 0, k = 0; i < fpdim; ++i) { 1307 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1308 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1309 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1310 } 1311 1312 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1313 PetscObject obj; 1314 PetscClassId id; 1315 PetscReal *B; 1316 PetscInt NcJ = 0, cpdim = 0, j; 1317 1318 ierr = PetscDSGetDiscretization(prob, fieldJ, &obj);CHKERRQ(ierr); 1319 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1320 if (id == PETSCFE_CLASSID) { 1321 PetscFE fe = (PetscFE) obj; 1322 1323 /* Evaluate basis at points */ 1324 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 1325 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1326 /* For now, fields only interpolate themselves */ 1327 if (fieldI == fieldJ) { 1328 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); 1329 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1330 for (i = 0, k = 0; i < fpdim; ++i) { 1331 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1332 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1333 for (p = 0; p < Np; ++p, ++k) { 1334 for (j = 0; j < cpdim; ++j) { 1335 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]; 1336 } 1337 } 1338 } 1339 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1340 } 1341 } else if (id == PETSCFV_CLASSID) { 1342 PetscFV fv = (PetscFV) obj; 1343 1344 /* Evaluate constant function at points */ 1345 ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr); 1346 cpdim = 1; 1347 /* For now, fields only interpolate themselves */ 1348 if (fieldI == fieldJ) { 1349 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); 1350 for (i = 0, k = 0; i < fpdim; ++i) { 1351 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1352 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1353 for (p = 0; p < Np; ++p, ++k) { 1354 for (j = 0; j < cpdim; ++j) { 1355 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p]; 1356 } 1357 } 1358 } 1359 } 1360 } 1361 offsetJ += cpdim*NcJ; 1362 } 1363 offsetI += fpdim*Nc; 1364 ierr = PetscFree(points);CHKERRQ(ierr); 1365 } 1366 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1367 /* Preallocate matrix */ 1368 { 1369 PetscHashJK ht; 1370 PetscLayout rLayout; 1371 PetscInt *dnz, *onz, *cellCIndices, *cellFIndices; 1372 PetscInt locRows, rStart, rEnd, cell, r; 1373 1374 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1375 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1376 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1377 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1378 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1379 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1380 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1381 ierr = PetscCalloc4(locRows,&dnz,locRows,&onz,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1382 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1383 for (cell = cStart; cell < cEnd; ++cell) { 1384 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1385 for (r = 0; r < rTotDim; ++r) { 1386 PetscHashJKKey key; 1387 PetscHashJKIter missing, iter; 1388 1389 key.j = cellFIndices[r]; 1390 if (key.j < 0) continue; 1391 for (c = 0; c < cTotDim; ++c) { 1392 key.k = cellCIndices[c]; 1393 if (key.k < 0) continue; 1394 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1395 if (missing) { 1396 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1397 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1398 else ++onz[key.j-rStart]; 1399 } 1400 } 1401 } 1402 } 1403 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1404 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1405 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1406 ierr = PetscFree4(dnz,onz,cellCIndices,cellFIndices);CHKERRQ(ierr); 1407 } 1408 /* Fill matrix */ 1409 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1410 for (c = cStart; c < cEnd; ++c) { 1411 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1412 } 1413 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1414 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1415 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1416 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1417 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1418 if (mesh->printFEM) { 1419 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1420 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1421 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1422 } 1423 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1424 PetscFunctionReturn(0); 1425 } 1426 1427 #undef __FUNCT__ 1428 #define __FUNCT__ "DMPlexComputeInterpolatorGeneral" 1429 /*@ 1430 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 1431 1432 Input Parameters: 1433 + dmf - The fine mesh 1434 . dmc - The coarse mesh 1435 - user - The user context 1436 1437 Output Parameter: 1438 . In - The interpolation matrix 1439 1440 Level: developer 1441 1442 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1443 @*/ 1444 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 1445 { 1446 PetscDS prob; 1447 PetscSection fsection, csection, globalFSection, globalCSection; 1448 PetscHashJK ht; 1449 PetscLayout rLayout; 1450 PetscInt *dnz, *onz; 1451 PetscInt locRows, rStart, rEnd; 1452 PetscReal *x, *v0, *J, *invJ, detJ; 1453 PetscReal *v0c, *Jc, *invJc, detJc; 1454 PetscScalar *elemMat; 1455 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1456 PetscErrorCode ierr; 1457 1458 PetscFunctionBegin; 1459 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1460 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1461 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1462 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1463 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1464 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1465 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1466 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1467 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1468 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1469 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1470 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1471 ierr = PetscMalloc1(totDim*totDim, &elemMat);CHKERRQ(ierr); 1472 1473 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1474 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1475 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1476 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1477 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1478 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1479 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1480 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1481 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1482 for (field = 0; field < Nf; ++field) { 1483 PetscObject obj; 1484 PetscClassId id; 1485 PetscDualSpace Q; 1486 PetscQuadrature f; 1487 const PetscReal *qpoints, *qweights; 1488 PetscInt npoints = 0, Nc, Np, fpdim, i, p, d; 1489 1490 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1491 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1492 if (id == PETSCFE_CLASSID) { 1493 PetscFE fe = (PetscFE) obj; 1494 1495 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1496 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1497 } else if (id == PETSCFV_CLASSID) { 1498 PetscFV fv = (PetscFV) obj; 1499 1500 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1501 Nc = 1; 1502 } 1503 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1504 /* For each fine grid cell */ 1505 for (cell = cStart; cell < cEnd; ++cell) { 1506 PetscInt *findices, *cindices; 1507 PetscInt numFIndices, numCIndices; 1508 1509 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1510 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1511 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1512 for (i = 0; i < fpdim; ++i) { 1513 Vec pointVec; 1514 PetscScalar *pV; 1515 IS coarseCellIS; 1516 const PetscInt *coarseCells; 1517 PetscInt numCoarseCells, q, r, c; 1518 1519 /* Get points from the dual basis functional quadrature */ 1520 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1521 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1522 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1523 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1524 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1525 for (q = 0; q < Np; ++q) { 1526 /* Transform point to real space */ 1527 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1528 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1529 } 1530 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1531 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1532 ierr = DMLocatePoints(dmc, pointVec, &coarseCellIS);CHKERRQ(ierr); 1533 ierr = ISViewFromOptions(coarseCellIS, NULL, "-interp_is_view");CHKERRQ(ierr); 1534 /* Update preallocation info */ 1535 ierr = ISGetLocalSize(coarseCellIS, &numCoarseCells);CHKERRQ(ierr); 1536 ierr = ISGetIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1537 for (r = 0; r < Nc; ++r) { 1538 PetscHashJKKey key; 1539 PetscHashJKIter missing, iter; 1540 1541 key.j = findices[i*Nc+r]; 1542 if (key.j < 0) continue; 1543 /* Get indices for coarse elements */ 1544 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1545 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1546 for (c = 0; c < numCIndices; ++c) { 1547 key.k = cindices[c]; 1548 if (key.k < 0) continue; 1549 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1550 if (missing) { 1551 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1552 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1553 else ++onz[key.j-rStart]; 1554 } 1555 } 1556 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1557 } 1558 } 1559 ierr = ISRestoreIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1560 ierr = ISDestroy(&coarseCellIS);CHKERRQ(ierr); 1561 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1562 } 1563 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1564 } 1565 } 1566 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1567 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1568 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1569 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1570 for (field = 0; field < Nf; ++field) { 1571 PetscObject obj; 1572 PetscClassId id; 1573 PetscDualSpace Q; 1574 PetscQuadrature f; 1575 const PetscReal *qpoints, *qweights; 1576 PetscInt npoints = 0, Nc, Np, fpdim, i, p, d; 1577 1578 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1579 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1580 if (id == PETSCFE_CLASSID) { 1581 PetscFE fe = (PetscFE) obj; 1582 1583 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1584 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1585 } else if (id == PETSCFV_CLASSID) { 1586 PetscFV fv = (PetscFV) obj; 1587 1588 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1589 Nc = 1; 1590 } 1591 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1592 /* For each fine grid cell */ 1593 for (cell = cStart; cell < cEnd; ++cell) { 1594 PetscInt *findices, *cindices; 1595 PetscInt numFIndices, numCIndices; 1596 1597 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1598 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1599 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1600 for (i = 0; i < fpdim; ++i) { 1601 Vec pointVec; 1602 PetscScalar *pV; 1603 IS coarseCellIS; 1604 const PetscInt *coarseCells; 1605 PetscInt numCoarseCells, cpdim, q, r, c, j; 1606 1607 /* Get points from the dual basis functional quadrature */ 1608 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1609 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, &qweights);CHKERRQ(ierr); 1610 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1611 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1612 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1613 for (q = 0; q < Np; ++q) { 1614 /* Transform point to real space */ 1615 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1616 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1617 } 1618 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1619 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1620 ierr = DMLocatePoints(dmc, pointVec, &coarseCellIS);CHKERRQ(ierr); 1621 /* Update preallocation info */ 1622 ierr = ISGetLocalSize(coarseCellIS, &numCoarseCells);CHKERRQ(ierr); 1623 ierr = ISGetIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1624 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1625 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1626 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1627 /* Transform points from real space to coarse reference space */ 1628 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell], NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 1629 CoordinatesRealToRef(dim, dim, v0c, invJc, &pV[ccell*dim], x); 1630 1631 if (id == PETSCFE_CLASSID) { 1632 PetscFE fe = (PetscFE) obj; 1633 PetscReal *B; 1634 1635 /* Evaluate coarse basis on contained point */ 1636 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1637 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 1638 /* Get elemMat entries by multiplying by weight */ 1639 for (j = 0; j < cpdim; ++j) { 1640 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = B[j*Nc + c]*qweights[ccell]; 1641 } 1642 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1643 } else { 1644 cpdim = 1; 1645 for (j = 0; j < cpdim; ++j) { 1646 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = 1.0*qweights[ccell]; 1647 } 1648 } 1649 /* Update interpolator */ 1650 ierr = MatSetValues(In, Nc, &findices[i*Nc], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1651 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);CHKERRQ(ierr);CHKERRQ(ierr); 1652 } 1653 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1654 ierr = ISRestoreIndices(coarseCellIS, &coarseCells);CHKERRQ(ierr); 1655 ierr = ISDestroy(&coarseCellIS);CHKERRQ(ierr); 1656 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1657 } 1658 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);CHKERRQ(ierr);CHKERRQ(ierr); 1659 } 1660 } 1661 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 1662 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 1663 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1664 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1665 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1666 PetscFunctionReturn(0); 1667 } 1668 1669 #undef __FUNCT__ 1670 #define __FUNCT__ "DMPlexComputeInjectorFEM" 1671 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 1672 { 1673 PetscDS prob; 1674 PetscFE *feRef; 1675 PetscFV *fvRef; 1676 Vec fv, cv; 1677 IS fis, cis; 1678 PetscSection fsection, fglobalSection, csection, cglobalSection; 1679 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 1680 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, offsetC, offsetF, m; 1681 PetscErrorCode ierr; 1682 1683 PetscFunctionBegin; 1684 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1685 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1686 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1687 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1688 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1689 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1690 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1691 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1692 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1693 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1694 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1695 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1696 for (f = 0; f < Nf; ++f) { 1697 PetscObject obj; 1698 PetscClassId id; 1699 PetscInt fNb = 0, Nc = 0; 1700 1701 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1702 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1703 if (id == PETSCFE_CLASSID) { 1704 PetscFE fe = (PetscFE) obj; 1705 1706 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1707 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 1708 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1709 } else if (id == PETSCFV_CLASSID) { 1710 PetscFV fv = (PetscFV) obj; 1711 PetscDualSpace Q; 1712 1713 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1714 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1715 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 1716 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1717 } 1718 fTotDim += fNb*Nc; 1719 } 1720 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1721 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 1722 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 1723 PetscFE feC; 1724 PetscFV fvC; 1725 PetscDualSpace QF, QC; 1726 PetscInt NcF, NcC, fpdim, cpdim; 1727 1728 if (feRef[field]) { 1729 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 1730 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 1731 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 1732 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 1733 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1734 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 1735 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1736 } else { 1737 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 1738 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 1739 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 1740 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 1741 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1742 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 1743 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1744 } 1745 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); 1746 for (c = 0; c < cpdim; ++c) { 1747 PetscQuadrature cfunc; 1748 const PetscReal *cqpoints; 1749 PetscInt NpC; 1750 PetscBool found = PETSC_FALSE; 1751 1752 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 1753 ierr = PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 1754 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 1755 for (f = 0; f < fpdim; ++f) { 1756 PetscQuadrature ffunc; 1757 const PetscReal *fqpoints; 1758 PetscReal sum = 0.0; 1759 PetscInt NpF, comp; 1760 1761 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 1762 ierr = PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 1763 if (NpC != NpF) continue; 1764 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 1765 if (sum > 1.0e-9) continue; 1766 for (comp = 0; comp < NcC; ++comp) { 1767 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 1768 } 1769 found = PETSC_TRUE; 1770 break; 1771 } 1772 if (!found) { 1773 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 1774 if (fvRef[field]) { 1775 PetscInt comp; 1776 for (comp = 0; comp < NcC; ++comp) { 1777 cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp; 1778 } 1779 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 1780 } 1781 } 1782 offsetC += cpdim*NcC; 1783 offsetF += fpdim*NcF; 1784 } 1785 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 1786 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1787 1788 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 1789 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 1790 ierr = VecGetOwnershipRange(cv, &startC, NULL);CHKERRQ(ierr); 1791 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 1792 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 1793 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 1794 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 1795 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 1796 for (c = cStart; c < cEnd; ++c) { 1797 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 1798 for (d = 0; d < cTotDim; ++d) { 1799 if (cellCIndices[d] < 0) continue; 1800 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]]); 1801 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 1802 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 1803 } 1804 } 1805 ierr = PetscFree(cmap);CHKERRQ(ierr); 1806 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 1807 1808 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 1809 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 1810 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 1811 ierr = ISDestroy(&cis);CHKERRQ(ierr); 1812 ierr = ISDestroy(&fis);CHKERRQ(ierr); 1813 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 1814 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 1815 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1816 PetscFunctionReturn(0); 1817 } 1818