1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petscsf.h> 3 #include <petscsnes.h> 4 5 #include <petsc/private/hashsetij.h> 6 #include <petsc/private/petscfeimpl.h> 7 #include <petsc/private/petscfvimpl.h> 8 9 static PetscErrorCode DMPlexConvertPlex(DM dm, DM *plex, PetscBool copy) 10 { 11 PetscBool isPlex; 12 PetscErrorCode ierr; 13 14 PetscFunctionBegin; 15 ierr = PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);CHKERRQ(ierr); 16 if (isPlex) { 17 *plex = dm; 18 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 19 } else { 20 ierr = PetscObjectQuery((PetscObject) dm, "dm_plex", (PetscObject *) plex);CHKERRQ(ierr); 21 if (!*plex) { 22 ierr = DMConvert(dm,DMPLEX,plex);CHKERRQ(ierr); 23 ierr = PetscObjectCompose((PetscObject) dm, "dm_plex", (PetscObject) *plex);CHKERRQ(ierr); 24 if (copy) { 25 const char *comps[3] = {"A", "dmAux"}; 26 PetscObject obj; 27 PetscInt i; 28 29 { 30 /* Run the subdomain hook (this will copy the DMSNES/DMTS) */ 31 DMSubDomainHookLink link; 32 for (link = dm->subdomainhook; link; link = link->next) { 33 if (link->ddhook) {ierr = (*link->ddhook)(dm, *plex, link->ctx);CHKERRQ(ierr);} 34 } 35 } 36 for (i = 0; i < 3; i++) { 37 ierr = PetscObjectQuery((PetscObject) dm, comps[i], &obj);CHKERRQ(ierr); 38 ierr = PetscObjectCompose((PetscObject) *plex, comps[i], obj);CHKERRQ(ierr); 39 } 40 } 41 } else { 42 ierr = PetscObjectReference((PetscObject) *plex);CHKERRQ(ierr); 43 } 44 } 45 PetscFunctionReturn(0); 46 } 47 48 static PetscErrorCode PetscContainerUserDestroy_PetscFEGeom (void *ctx) 49 { 50 PetscFEGeom *geom = (PetscFEGeom *) ctx; 51 PetscErrorCode ierr; 52 53 PetscFunctionBegin; 54 ierr = PetscFEGeomDestroy(&geom);CHKERRQ(ierr); 55 PetscFunctionReturn(0); 56 } 57 58 static PetscErrorCode DMPlexGetFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom) 59 { 60 char composeStr[33] = {0}; 61 PetscObjectId id; 62 PetscContainer container; 63 PetscErrorCode ierr; 64 65 PetscFunctionBegin; 66 ierr = PetscObjectGetId((PetscObject)quad,&id);CHKERRQ(ierr); 67 ierr = PetscSNPrintf(composeStr, 32, "DMPlexGetFEGeom_%x\n", id);CHKERRQ(ierr); 68 ierr = PetscObjectQuery((PetscObject) pointIS, composeStr, (PetscObject *) &container);CHKERRQ(ierr); 69 if (container) { 70 ierr = PetscContainerGetPointer(container, (void **) geom);CHKERRQ(ierr); 71 } else { 72 ierr = DMFieldCreateFEGeom(coordField, pointIS, quad, faceData, geom);CHKERRQ(ierr); 73 ierr = PetscContainerCreate(PETSC_COMM_SELF,&container);CHKERRQ(ierr); 74 ierr = PetscContainerSetPointer(container, (void *) *geom);CHKERRQ(ierr); 75 ierr = PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom);CHKERRQ(ierr); 76 ierr = PetscObjectCompose((PetscObject) pointIS, composeStr, (PetscObject) container);CHKERRQ(ierr); 77 ierr = PetscContainerDestroy(&container);CHKERRQ(ierr); 78 } 79 PetscFunctionReturn(0); 80 } 81 82 static PetscErrorCode DMPlexRestoreFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom) 83 { 84 PetscFunctionBegin; 85 *geom = NULL; 86 PetscFunctionReturn(0); 87 } 88 89 /*@ 90 DMPlexGetScale - Get the scale for the specified fundamental unit 91 92 Not collective 93 94 Input Arguments: 95 + dm - the DM 96 - unit - The SI unit 97 98 Output Argument: 99 . scale - The value used to scale all quantities with this unit 100 101 Level: advanced 102 103 .seealso: DMPlexSetScale(), PetscUnit 104 @*/ 105 PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale) 106 { 107 DM_Plex *mesh = (DM_Plex*) dm->data; 108 109 PetscFunctionBegin; 110 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 111 PetscValidPointer(scale, 3); 112 *scale = mesh->scale[unit]; 113 PetscFunctionReturn(0); 114 } 115 116 /*@ 117 DMPlexSetScale - Set the scale for the specified fundamental unit 118 119 Not collective 120 121 Input Arguments: 122 + dm - the DM 123 . unit - The SI unit 124 - scale - The value used to scale all quantities with this unit 125 126 Level: advanced 127 128 .seealso: DMPlexGetScale(), PetscUnit 129 @*/ 130 PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale) 131 { 132 DM_Plex *mesh = (DM_Plex*) dm->data; 133 134 PetscFunctionBegin; 135 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 136 mesh->scale[unit] = scale; 137 PetscFunctionReturn(0); 138 } 139 140 static PetscErrorCode DMPlexProjectRigidBody_Private(PetscInt dim, PetscReal t, const PetscReal X[], PetscInt Nf, PetscScalar *mode, void *ctx) 141 { 142 const PetscInt eps[3][3][3] = {{{0, 0, 0}, {0, 0, 1}, {0, -1, 0}}, {{0, 0, -1}, {0, 0, 0}, {1, 0, 0}}, {{0, 1, 0}, {-1, 0, 0}, {0, 0, 0}}}; 143 PetscInt *ctxInt = (PetscInt *) ctx; 144 PetscInt dim2 = ctxInt[0]; 145 PetscInt d = ctxInt[1]; 146 PetscInt i, j, k = dim > 2 ? d - dim : d; 147 148 PetscFunctionBegin; 149 if (dim != dim2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input dimension %d does not match context dimension %d", dim, dim2); 150 for (i = 0; i < dim; i++) mode[i] = 0.; 151 if (d < dim) { 152 mode[d] = 1.; /* Translation along axis d */ 153 } else { 154 for (i = 0; i < dim; i++) { 155 for (j = 0; j < dim; j++) { 156 mode[j] += eps[i][j][k]*X[i]; /* Rotation about axis d */ 157 } 158 } 159 } 160 PetscFunctionReturn(0); 161 } 162 163 /*@ 164 DMPlexCreateRigidBody - For the default global section, create rigid body modes by function space interpolation 165 166 Collective on DM 167 168 Input Arguments: 169 . dm - the DM 170 171 Output Argument: 172 . sp - the null space 173 174 Note: This is necessary to provide a suitable coarse space for algebraic multigrid 175 176 Level: advanced 177 178 .seealso: MatNullSpaceCreate(), PCGAMG 179 @*/ 180 PetscErrorCode DMPlexCreateRigidBody(DM dm, MatNullSpace *sp) 181 { 182 MPI_Comm comm; 183 Vec mode[6]; 184 PetscSection section, globalSection; 185 PetscInt dim, dimEmbed, n, m, mmin, d, i, j; 186 PetscErrorCode ierr; 187 188 PetscFunctionBegin; 189 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 190 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 191 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 192 if (dim == 1) { 193 ierr = MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);CHKERRQ(ierr); 194 PetscFunctionReturn(0); 195 } 196 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 197 ierr = DMGetGlobalSection(dm, &globalSection);CHKERRQ(ierr); 198 ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr); 199 m = (dim*(dim+1))/2; 200 ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr); 201 ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr); 202 ierr = VecSetUp(mode[0]);CHKERRQ(ierr); 203 ierr = VecGetSize(mode[0], &n);CHKERRQ(ierr); 204 mmin = PetscMin(m, n); 205 for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);} 206 for (d = 0; d < m; d++) { 207 PetscInt ctx[2]; 208 PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private; 209 void *voidctx = (void *) (&ctx[0]); 210 211 ctx[0] = dimEmbed; 212 ctx[1] = d; 213 ierr = DMProjectFunction(dm, 0.0, &func, &voidctx, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 214 } 215 for (i = 0; i < PetscMin(dim, mmin); ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);} 216 /* Orthonormalize system */ 217 for (i = dim; i < mmin; ++i) { 218 PetscScalar dots[6]; 219 220 ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr); 221 for (j = 0; j < i; ++j) dots[j] *= -1.0; 222 ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr); 223 ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr); 224 } 225 ierr = MatNullSpaceCreate(comm, PETSC_FALSE, mmin, mode, sp);CHKERRQ(ierr); 226 for (i = 0; i < m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);} 227 PetscFunctionReturn(0); 228 } 229 230 /*@ 231 DMPlexCreateRigidBodies - For the default global section, create rigid body modes by function space interpolation 232 233 Collective on DM 234 235 Input Arguments: 236 + dm - the DM 237 . nb - The number of bodies 238 . label - The DMLabel marking each domain 239 . nids - The number of ids per body 240 - ids - An array of the label ids in sequence for each domain 241 242 Output Argument: 243 . sp - the null space 244 245 Note: This is necessary to provide a suitable coarse space for algebraic multigrid 246 247 Level: advanced 248 249 .seealso: MatNullSpaceCreate() 250 @*/ 251 PetscErrorCode DMPlexCreateRigidBodies(DM dm, PetscInt nb, DMLabel label, const PetscInt nids[], const PetscInt ids[], MatNullSpace *sp) 252 { 253 MPI_Comm comm; 254 PetscSection section, globalSection; 255 Vec *mode; 256 PetscScalar *dots; 257 PetscInt dim, dimEmbed, n, m, b, d, i, j, off; 258 PetscErrorCode ierr; 259 260 PetscFunctionBegin; 261 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 262 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 263 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 264 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 265 ierr = DMGetGlobalSection(dm, &globalSection);CHKERRQ(ierr); 266 ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr); 267 m = nb * (dim*(dim+1))/2; 268 ierr = PetscMalloc2(m, &mode, m, &dots);CHKERRQ(ierr); 269 ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr); 270 ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr); 271 ierr = VecSetUp(mode[0]);CHKERRQ(ierr); 272 for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);} 273 for (b = 0, off = 0; b < nb; ++b) { 274 for (d = 0; d < m/nb; ++d) { 275 PetscInt ctx[2]; 276 PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private; 277 void *voidctx = (void *) (&ctx[0]); 278 279 ctx[0] = dimEmbed; 280 ctx[1] = d; 281 ierr = DMProjectFunctionLabel(dm, 0.0, label, nids[b], &ids[off], 0, NULL, &func, &voidctx, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 282 off += nids[b]; 283 } 284 } 285 for (i = 0; i < dim; ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);} 286 /* Orthonormalize system */ 287 for (i = 0; i < m; ++i) { 288 ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr); 289 for (j = 0; j < i; ++j) dots[j] *= -1.0; 290 ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr); 291 ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr); 292 } 293 ierr = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);CHKERRQ(ierr); 294 for (i = 0; i< m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);} 295 ierr = PetscFree2(mode, dots);CHKERRQ(ierr); 296 PetscFunctionReturn(0); 297 } 298 299 /*@ 300 DMPlexSetMaxProjectionHeight - In DMPlexProjectXXXLocal() functions, the projected values of a basis function's dofs 301 are computed by associating the basis function with one of the mesh points in its transitively-closed support, and 302 evaluating the dual space basis of that point. A basis function is associated with the point in its 303 transitively-closed support whose mesh height is highest (w.r.t. DAG height), but not greater than the maximum 304 projection height, which is set with this function. By default, the maximum projection height is zero, which means 305 that only mesh cells are used to project basis functions. A height of one, for example, evaluates a cell-interior 306 basis functions using its cells dual space basis, but all other basis functions with the dual space basis of a face. 307 308 Input Parameters: 309 + dm - the DMPlex object 310 - height - the maximum projection height >= 0 311 312 Level: advanced 313 314 .seealso: DMPlexGetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal() 315 @*/ 316 PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height) 317 { 318 DM_Plex *plex = (DM_Plex *) dm->data; 319 320 PetscFunctionBegin; 321 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 322 plex->maxProjectionHeight = height; 323 PetscFunctionReturn(0); 324 } 325 326 /*@ 327 DMPlexGetMaxProjectionHeight - Get the maximum height (w.r.t. DAG) of mesh points used to evaluate dual bases in 328 DMPlexProjectXXXLocal() functions. 329 330 Input Parameters: 331 . dm - the DMPlex object 332 333 Output Parameters: 334 . height - the maximum projection height 335 336 Level: intermediate 337 338 .seealso: DMPlexSetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal() 339 @*/ 340 PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height) 341 { 342 DM_Plex *plex = (DM_Plex *) dm->data; 343 344 PetscFunctionBegin; 345 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 346 *height = plex->maxProjectionHeight; 347 PetscFunctionReturn(0); 348 } 349 350 typedef struct { 351 PetscReal alpha; /* The first Euler angle, and in 2D the only one */ 352 PetscReal beta; /* The second Euler angle */ 353 PetscReal gamma; /* The third Euler angle */ 354 PetscInt dim; /* The dimension of R */ 355 PetscScalar *R; /* The rotation matrix, transforming a vector in the local basis to the global basis */ 356 PetscScalar *RT; /* The transposed rotation matrix, transforming a vector in the global basis to the local basis */ 357 } RotCtx; 358 359 /* 360 Note: Following https://en.wikipedia.org/wiki/Euler_angles, we will specify Euler angles by extrinsic rotations, meaning that 361 we rotate with respect to a fixed initial coordinate system, the local basis (x-y-z). The global basis (X-Y-Z) is reached as follows: 362 $ The XYZ system rotates about the z axis by alpha. The X axis is now at angle alpha with respect to the x axis. 363 $ The XYZ system rotates again about the x axis by beta. The Z axis is now at angle beta with respect to the z axis. 364 $ The XYZ system rotates a third time about the z axis by gamma. 365 */ 366 static PetscErrorCode DMPlexBasisTransformSetUp_Rotation_Internal(DM dm, void *ctx) 367 { 368 RotCtx *rc = (RotCtx *) ctx; 369 PetscInt dim = rc->dim; 370 PetscReal c1, s1, c2, s2, c3, s3; 371 PetscErrorCode ierr; 372 373 PetscFunctionBegin; 374 ierr = PetscMalloc2(PetscSqr(dim), &rc->R, PetscSqr(dim), &rc->RT);CHKERRQ(ierr); 375 switch (dim) { 376 case 2: 377 c1 = PetscCosReal(rc->alpha);s1 = PetscSinReal(rc->alpha); 378 rc->R[0] = c1;rc->R[1] = s1; 379 rc->R[2] = -s1;rc->R[3] = c1; 380 ierr = PetscMemcpy(rc->RT, rc->R, PetscSqr(dim) * sizeof(PetscScalar));CHKERRQ(ierr); 381 DMPlex_Transpose2D_Internal(rc->RT);break; 382 break; 383 case 3: 384 c1 = PetscCosReal(rc->alpha);s1 = PetscSinReal(rc->alpha); 385 c2 = PetscCosReal(rc->beta); s2 = PetscSinReal(rc->beta); 386 c3 = PetscCosReal(rc->gamma);s3 = PetscSinReal(rc->gamma); 387 rc->R[0] = c1*c3 - c2*s1*s3;rc->R[1] = c3*s1 + c1*c2*s3;rc->R[2] = s2*s3; 388 rc->R[3] = -c1*s3 - c2*c3*s1;rc->R[4] = c1*c2*c3 - s1*s3; rc->R[5] = c3*s2; 389 rc->R[6] = s1*s2; rc->R[7] = -c1*s2; rc->R[8] = c2; 390 ierr = PetscMemcpy(rc->RT, rc->R, PetscSqr(dim) * sizeof(PetscScalar));CHKERRQ(ierr); 391 DMPlex_Transpose3D_Internal(rc->RT);break; 392 break; 393 default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not supported", dim); 394 } 395 PetscFunctionReturn(0); 396 } 397 398 static PetscErrorCode DMPlexBasisTransformDestroy_Rotation_Internal(DM dm, void *ctx) 399 { 400 RotCtx *rc = (RotCtx *) ctx; 401 PetscErrorCode ierr; 402 403 PetscFunctionBegin; 404 ierr = PetscFree2(rc->R, rc->RT);CHKERRQ(ierr); 405 ierr = PetscFree(rc);CHKERRQ(ierr); 406 PetscFunctionReturn(0); 407 } 408 409 static PetscErrorCode DMPlexBasisTransformGetMatrix_Rotation_Internal(DM dm, const PetscReal x[], PetscBool l2g, const PetscScalar **A, void *ctx) 410 { 411 RotCtx *rc = (RotCtx *) ctx; 412 413 PetscFunctionBeginHot; 414 PetscValidPointer(ctx, 5); 415 if (l2g) {*A = rc->R;} 416 else {*A = rc->RT;} 417 PetscFunctionReturn(0); 418 } 419 420 PetscErrorCode DMPlexBasisTransformApplyReal_Internal(DM dm, const PetscReal x[], PetscBool l2g, PetscInt dim, const PetscReal *y, PetscReal *z, void *ctx) 421 { 422 PetscErrorCode ierr; 423 424 PetscFunctionBegin; 425 #if defined(PETSC_USE_COMPLEX) 426 switch (dim) { 427 case 2: 428 { 429 PetscScalar yt[2], zt[2]; 430 431 yt[0] = y[0]; yt[1] = y[1]; 432 ierr = DMPlexBasisTransformApply_Internal(dm, x, l2g, dim, yt, zt, ctx);CHKERRQ(ierr); 433 z[0] = PetscRealPart(zt[0]); z[1] = PetscRealPart(zt[1]); 434 } 435 break; 436 case 3: 437 { 438 PetscScalar yt[3], zt[3]; 439 440 yt[0] = y[0]; yt[1] = y[1]; yt[2] = y[2]; 441 ierr = DMPlexBasisTransformApply_Internal(dm, x, l2g, dim, yt, zt, ctx);CHKERRQ(ierr); 442 z[0] = PetscRealPart(zt[0]); z[1] = PetscRealPart(zt[1]); z[2] = PetscRealPart(zt[2]); 443 } 444 break; 445 } 446 #else 447 ierr = DMPlexBasisTransformApply_Internal(dm, x, l2g, dim, y, z, ctx);CHKERRQ(ierr); 448 #endif 449 PetscFunctionReturn(0); 450 } 451 452 PetscErrorCode DMPlexBasisTransformApply_Internal(DM dm, const PetscReal x[], PetscBool l2g, PetscInt dim, const PetscScalar *y, PetscScalar *z, void *ctx) 453 { 454 const PetscScalar *A; 455 PetscErrorCode ierr; 456 457 PetscFunctionBeginHot; 458 ierr = (*dm->transformGetMatrix)(dm, x, l2g, &A, ctx);CHKERRQ(ierr); 459 switch (dim) { 460 case 2: DMPlex_Mult2D_Internal(A, 1, y, z);break; 461 case 3: DMPlex_Mult3D_Internal(A, 1, y, z);break; 462 } 463 PetscFunctionReturn(0); 464 } 465 466 static PetscErrorCode DMPlexBasisTransformField_Internal(DM dm, DM tdm, Vec tv, PetscInt p, PetscInt f, PetscBool l2g, PetscScalar *a) 467 { 468 PetscSection ts; 469 const PetscScalar *ta, *tva; 470 PetscInt dof; 471 PetscErrorCode ierr; 472 473 PetscFunctionBeginHot; 474 ierr = DMGetSection(tdm, &ts);CHKERRQ(ierr); 475 ierr = PetscSectionGetFieldDof(ts, p, f, &dof);CHKERRQ(ierr); 476 ierr = VecGetArrayRead(tv, &ta);CHKERRQ(ierr); 477 ierr = DMPlexPointLocalFieldRead(tdm, p, f, ta, (void *) &tva);CHKERRQ(ierr); 478 if (l2g) { 479 switch (dof) { 480 case 4: DMPlex_Mult2D_Internal(tva, 1, a, a);break; 481 case 9: DMPlex_Mult3D_Internal(tva, 1, a, a);break; 482 } 483 } else { 484 switch (dof) { 485 case 4: DMPlex_MultTranspose2D_Internal(tva, 1, a, a);break; 486 case 9: DMPlex_MultTranspose3D_Internal(tva, 1, a, a);break; 487 } 488 } 489 ierr = VecRestoreArrayRead(tv, &ta);CHKERRQ(ierr); 490 PetscFunctionReturn(0); 491 } 492 493 static PetscErrorCode DMPlexBasisTransformFieldTensor_Internal(DM dm, DM tdm, Vec tv, PetscInt pf, PetscInt f, PetscInt pg, PetscInt g, PetscBool l2g, PetscInt lda, PetscScalar *a) 494 { 495 PetscSection s, ts; 496 const PetscScalar *ta, *tvaf, *tvag; 497 PetscInt fdof, gdof, fpdof, gpdof; 498 PetscErrorCode ierr; 499 500 PetscFunctionBeginHot; 501 ierr = DMGetSection(dm, &s);CHKERRQ(ierr); 502 ierr = DMGetSection(tdm, &ts);CHKERRQ(ierr); 503 ierr = PetscSectionGetFieldDof(s, pf, f, &fpdof);CHKERRQ(ierr); 504 ierr = PetscSectionGetFieldDof(s, pg, g, &gpdof);CHKERRQ(ierr); 505 ierr = PetscSectionGetFieldDof(ts, pf, f, &fdof);CHKERRQ(ierr); 506 ierr = PetscSectionGetFieldDof(ts, pg, g, &gdof);CHKERRQ(ierr); 507 ierr = VecGetArrayRead(tv, &ta);CHKERRQ(ierr); 508 ierr = DMPlexPointLocalFieldRead(tdm, pf, f, ta, (void *) &tvaf);CHKERRQ(ierr); 509 ierr = DMPlexPointLocalFieldRead(tdm, pg, g, ta, (void *) &tvag);CHKERRQ(ierr); 510 if (l2g) { 511 switch (fdof) { 512 case 4: DMPlex_MatMult2D_Internal(tvaf, gpdof, lda, a, a);break; 513 case 9: DMPlex_MatMult3D_Internal(tvaf, gpdof, lda, a, a);break; 514 } 515 switch (gdof) { 516 case 4: DMPlex_MatMultTransposeLeft2D_Internal(tvag, fpdof, lda, a, a);break; 517 case 9: DMPlex_MatMultTransposeLeft3D_Internal(tvag, fpdof, lda, a, a);break; 518 } 519 } else { 520 switch (fdof) { 521 case 4: DMPlex_MatMultTranspose2D_Internal(tvaf, gpdof, lda, a, a);break; 522 case 9: DMPlex_MatMultTranspose3D_Internal(tvaf, gpdof, lda, a, a);break; 523 } 524 switch (gdof) { 525 case 4: DMPlex_MatMultLeft2D_Internal(tvag, fpdof, lda, a, a);break; 526 case 9: DMPlex_MatMultLeft3D_Internal(tvag, fpdof, lda, a, a);break; 527 } 528 } 529 ierr = VecRestoreArrayRead(tv, &ta);CHKERRQ(ierr); 530 PetscFunctionReturn(0); 531 } 532 533 PetscErrorCode DMPlexBasisTransformPoint_Internal(DM dm, DM tdm, Vec tv, PetscInt p, PetscBool fieldActive[], PetscBool l2g, PetscScalar *a) 534 { 535 PetscSection s; 536 PetscSection clSection; 537 IS clPoints; 538 const PetscInt *clp; 539 PetscInt *points = NULL; 540 PetscInt Nf, f, Np, cp, dof, d = 0; 541 PetscErrorCode ierr; 542 543 PetscFunctionBegin; 544 ierr = DMGetSection(dm, &s);CHKERRQ(ierr); 545 ierr = PetscSectionGetNumFields(s, &Nf);CHKERRQ(ierr); 546 ierr = DMPlexGetCompressedClosure(dm, s, p, &Np, &points, &clSection, &clPoints, &clp);CHKERRQ(ierr); 547 for (f = 0; f < Nf; ++f) { 548 for (cp = 0; cp < Np*2; cp += 2) { 549 ierr = PetscSectionGetFieldDof(s, points[cp], f, &dof);CHKERRQ(ierr); 550 if (!dof) continue; 551 if (fieldActive[f]) {ierr = DMPlexBasisTransformField_Internal(dm, tdm, tv, points[cp], f, l2g, &a[d]);CHKERRQ(ierr);} 552 d += dof; 553 } 554 } 555 ierr = DMPlexRestoreCompressedClosure(dm, s, p, &Np, &points, &clSection, &clPoints, &clp);CHKERRQ(ierr); 556 PetscFunctionReturn(0); 557 } 558 559 PetscErrorCode DMPlexBasisTransformPointTensor_Internal(DM dm, DM tdm, Vec tv, PetscInt p, PetscBool l2g, PetscInt lda, PetscScalar *a) 560 { 561 PetscSection s; 562 PetscSection clSection; 563 IS clPoints; 564 const PetscInt *clp; 565 PetscInt *points = NULL; 566 PetscInt Nf, f, g, Np, cpf, cpg, fdof, gdof, r, c = 0; 567 PetscErrorCode ierr; 568 569 PetscFunctionBegin; 570 ierr = DMGetSection(dm, &s);CHKERRQ(ierr); 571 ierr = PetscSectionGetNumFields(s, &Nf);CHKERRQ(ierr); 572 ierr = DMPlexGetCompressedClosure(dm, s, p, &Np, &points, &clSection, &clPoints, &clp);CHKERRQ(ierr); 573 for (f = 0, r = 0; f < Nf; ++f) { 574 for (cpf = 0; cpf < Np*2; cpf += 2) { 575 ierr = PetscSectionGetFieldDof(s, points[cpf], f, &fdof);CHKERRQ(ierr); 576 for (g = 0, c = 0; g < Nf; ++g) { 577 for (cpg = 0; cpg < Np*2; cpg += 2) { 578 ierr = PetscSectionGetFieldDof(s, points[cpg], g, &gdof);CHKERRQ(ierr); 579 ierr = DMPlexBasisTransformFieldTensor_Internal(dm, tdm, tv, points[cpf], f, points[cpg], g, l2g, lda, &a[r*lda+c]);CHKERRQ(ierr); 580 c += gdof; 581 } 582 } 583 if (c != lda) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of columns %D should be %D", c, lda); 584 r += fdof; 585 } 586 } 587 if (r != lda) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of rows %D should be %D", c, lda); 588 ierr = DMPlexRestoreCompressedClosure(dm, s, p, &Np, &points, &clSection, &clPoints, &clp);CHKERRQ(ierr); 589 PetscFunctionReturn(0); 590 } 591 592 static PetscErrorCode DMPlexBasisTransform_Internal(DM dm, Vec lv, PetscBool l2g) 593 { 594 DM tdm; 595 Vec tv; 596 PetscSection ts, s; 597 const PetscScalar *ta; 598 PetscScalar *a, *va; 599 PetscInt pStart, pEnd, p, Nf, f; 600 PetscErrorCode ierr; 601 602 PetscFunctionBegin; 603 ierr = DMGetBasisTransformDM_Internal(dm, &tdm);CHKERRQ(ierr); 604 ierr = DMGetBasisTransformVec_Internal(dm, &tv);CHKERRQ(ierr); 605 ierr = DMGetSection(tdm, &ts);CHKERRQ(ierr); 606 ierr = DMGetSection(dm, &s);CHKERRQ(ierr); 607 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr); 608 ierr = PetscSectionGetNumFields(s, &Nf);CHKERRQ(ierr); 609 ierr = VecGetArray(lv, &a);CHKERRQ(ierr); 610 ierr = VecGetArrayRead(tv, &ta);CHKERRQ(ierr); 611 for (p = pStart; p < pEnd; ++p) { 612 for (f = 0; f < Nf; ++f) { 613 ierr = DMPlexPointLocalFieldRef(dm, p, f, a, (void *) &va);CHKERRQ(ierr); 614 ierr = DMPlexBasisTransformField_Internal(dm, tdm, tv, p, f, l2g, va);CHKERRQ(ierr); 615 } 616 } 617 ierr = VecRestoreArray(lv, &a);CHKERRQ(ierr); 618 ierr = VecRestoreArrayRead(tv, &ta);CHKERRQ(ierr); 619 PetscFunctionReturn(0); 620 } 621 622 /*@ 623 DMPlexGlobalToLocalBasis - Transform the values in the given local vector from the global basis to the local basis 624 625 Input Parameters: 626 + dm - The DM 627 - lv - A local vector with values in the global basis 628 629 Output Parameters: 630 . lv - A local vector with values in the local basis 631 632 Note: This method is only intended to be called inside DMGlobalToLocal(). It is unlikely that a user will have a local vector full of coefficients for the global basis unless they are reimplementing GlobalToLocal. 633 634 Level: developer 635 636 .seealso: DMPlexLocalToGlobalBasis(), DMGetSection() 637 @*/ 638 PetscErrorCode DMPlexGlobalToLocalBasis(DM dm, Vec lv) 639 { 640 PetscErrorCode ierr; 641 642 PetscFunctionBegin; 643 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 644 PetscValidHeaderSpecific(lv, VEC_CLASSID, 2); 645 ierr = DMPlexBasisTransform_Internal(dm, lv, PETSC_FALSE);CHKERRQ(ierr); 646 PetscFunctionReturn(0); 647 } 648 649 /*@ 650 DMPlexLocalToGlobalBasis - Transform the values in the given local vector from the local basis to the global basis 651 652 Input Parameters: 653 + dm - The DM 654 - lv - A local vector with values in the local basis 655 656 Output Parameters: 657 . lv - A local vector with values in the global basis 658 659 Note: This method is only intended to be called inside DMGlobalToLocal(). It is unlikely that a user would want a local vector full of coefficients for the global basis unless they are reimplementing GlobalToLocal. 660 661 Level: developer 662 663 .seealso: DMPlexGlobalToLocalBasis(), DMGetSection() 664 @*/ 665 PetscErrorCode DMPlexLocalToGlobalBasis(DM dm, Vec lv) 666 { 667 PetscErrorCode ierr; 668 669 PetscFunctionBegin; 670 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 671 PetscValidHeaderSpecific(lv, VEC_CLASSID, 2); 672 ierr = DMPlexBasisTransform_Internal(dm, lv, PETSC_TRUE);CHKERRQ(ierr); 673 PetscFunctionReturn(0); 674 } 675 676 /*@ 677 DMPlexCreateBasisRotation - Create an internal transformation from the global basis, used to specify boundary conditions 678 and global solutions, to a local basis, appropriate for discretization integrals and assembly. 679 680 Input Parameters: 681 + dm - The DM 682 . alpha - The first Euler angle, and in 2D the only one 683 . beta - The second Euler angle 684 . gamma - The third Euler angle 685 686 Note: Following https://en.wikipedia.org/wiki/Euler_angles, we will specify Euler angles by extrinsic rotations, meaning that 687 we rotate with respect to a fixed initial coordinate system, the local basis (x-y-z). The global basis (X-Y-Z) is reached as follows: 688 $ The XYZ system rotates about the z axis by alpha. The X axis is now at angle alpha with respect to the x axis. 689 $ The XYZ system rotates again about the x axis by beta. The Z axis is now at angle beta with respect to the z axis. 690 $ The XYZ system rotates a third time about the z axis by gamma. 691 692 Level: developer 693 694 .seealso: DMPlexGlobalToLocalBasis(), DMPlexLocalToGlobalBasis() 695 @*/ 696 PetscErrorCode DMPlexCreateBasisRotation(DM dm, PetscReal alpha, PetscReal beta, PetscReal gamma) 697 { 698 RotCtx *rc; 699 PetscInt cdim; 700 PetscErrorCode ierr; 701 702 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 703 ierr = PetscMalloc1(1, &rc);CHKERRQ(ierr); 704 dm->transformCtx = rc; 705 dm->transformSetUp = DMPlexBasisTransformSetUp_Rotation_Internal; 706 dm->transformDestroy = DMPlexBasisTransformDestroy_Rotation_Internal; 707 dm->transformGetMatrix = DMPlexBasisTransformGetMatrix_Rotation_Internal; 708 rc->dim = cdim; 709 rc->alpha = alpha; 710 rc->beta = beta; 711 rc->gamma = gamma; 712 ierr = (*dm->transformSetUp)(dm, dm->transformCtx);CHKERRQ(ierr); 713 ierr = DMConstructBasisTransform_Internal(dm);CHKERRQ(ierr); 714 PetscFunctionReturn(0); 715 } 716 717 /*@C 718 DMPlexInsertBoundaryValuesEssential - Insert boundary values into a local vector 719 720 Input Parameters: 721 + dm - The DM, with a PetscDS that matches the problem being constrained 722 . time - The time 723 . field - The field to constrain 724 . Nc - The number of constrained field components, or 0 for all components 725 . comps - An array of constrained component numbers, or NULL for all components 726 . label - The DMLabel defining constrained points 727 . numids - The number of DMLabel ids for constrained points 728 . ids - An array of ids for constrained points 729 . func - A pointwise function giving boundary values 730 - ctx - An optional user context for bcFunc 731 732 Output Parameter: 733 . locX - A local vector to receives the boundary values 734 735 Level: developer 736 737 .seealso: DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary() 738 @*/ 739 PetscErrorCode DMPlexInsertBoundaryValuesEssential(DM dm, PetscReal time, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void *ctx, Vec locX) 740 { 741 PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx); 742 void **ctxs; 743 PetscInt numFields; 744 PetscErrorCode ierr; 745 746 PetscFunctionBegin; 747 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 748 ierr = PetscCalloc2(numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 749 funcs[field] = func; 750 ctxs[field] = ctx; 751 ierr = DMProjectFunctionLabelLocal(dm, time, label, numids, ids, Nc, comps, funcs, ctxs, INSERT_BC_VALUES, locX);CHKERRQ(ierr); 752 ierr = PetscFree2(funcs,ctxs);CHKERRQ(ierr); 753 PetscFunctionReturn(0); 754 } 755 756 /*@C 757 DMPlexInsertBoundaryValuesEssentialField - Insert boundary values into a local vector 758 759 Input Parameters: 760 + dm - The DM, with a PetscDS that matches the problem being constrained 761 . time - The time 762 . locU - A local vector with the input solution values 763 . field - The field to constrain 764 . Nc - The number of constrained field components, or 0 for all components 765 . comps - An array of constrained component numbers, or NULL for all components 766 . label - The DMLabel defining constrained points 767 . numids - The number of DMLabel ids for constrained points 768 . ids - An array of ids for constrained points 769 . func - A pointwise function giving boundary values 770 - ctx - An optional user context for bcFunc 771 772 Output Parameter: 773 . locX - A local vector to receives the boundary values 774 775 Level: developer 776 777 .seealso: DMPlexInsertBoundaryValuesEssential(), DMAddBoundary() 778 @*/ 779 PetscErrorCode DMPlexInsertBoundaryValuesEssentialField(DM dm, PetscReal time, Vec locU, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[], 780 void (*func)(PetscInt, PetscInt, PetscInt, 781 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 782 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 783 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], 784 PetscScalar[]), 785 void *ctx, Vec locX) 786 { 787 void (**funcs)(PetscInt, PetscInt, PetscInt, 788 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 789 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 790 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]); 791 void **ctxs; 792 PetscInt numFields; 793 PetscErrorCode ierr; 794 795 PetscFunctionBegin; 796 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 797 ierr = PetscCalloc2(numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 798 funcs[field] = func; 799 ctxs[field] = ctx; 800 ierr = DMProjectFieldLabelLocal(dm, time, label, numids, ids, Nc, comps, locU, funcs, INSERT_BC_VALUES, locX);CHKERRQ(ierr); 801 ierr = PetscFree2(funcs,ctxs);CHKERRQ(ierr); 802 PetscFunctionReturn(0); 803 } 804 805 /*@C 806 DMPlexInsertBoundaryValuesRiemann - Insert boundary values into a local vector 807 808 Input Parameters: 809 + dm - The DM, with a PetscDS that matches the problem being constrained 810 . time - The time 811 . faceGeometry - A vector with the FVM face geometry information 812 . cellGeometry - A vector with the FVM cell geometry information 813 . Grad - A vector with the FVM cell gradient information 814 . field - The field to constrain 815 . Nc - The number of constrained field components, or 0 for all components 816 . comps - An array of constrained component numbers, or NULL for all components 817 . label - The DMLabel defining constrained points 818 . numids - The number of DMLabel ids for constrained points 819 . ids - An array of ids for constrained points 820 . func - A pointwise function giving boundary values 821 - ctx - An optional user context for bcFunc 822 823 Output Parameter: 824 . locX - A local vector to receives the boundary values 825 826 Note: This implementation currently ignores the numcomps/comps argument from DMAddBoundary() 827 828 Level: developer 829 830 .seealso: DMPlexInsertBoundaryValuesEssential(), DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary() 831 @*/ 832 PetscErrorCode DMPlexInsertBoundaryValuesRiemann(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[], 833 PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX) 834 { 835 PetscDS prob; 836 PetscSF sf; 837 DM dmFace, dmCell, dmGrad; 838 const PetscScalar *facegeom, *cellgeom = NULL, *grad; 839 const PetscInt *leaves; 840 PetscScalar *x, *fx; 841 PetscInt dim, nleaves, loc, fStart, fEnd, pdim, i; 842 PetscErrorCode ierr, ierru = 0; 843 844 PetscFunctionBegin; 845 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 846 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 847 nleaves = PetscMax(0, nleaves); 848 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 849 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 850 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 851 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 852 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 853 if (cellGeometry) { 854 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 855 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 856 } 857 if (Grad) { 858 PetscFV fv; 859 860 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);CHKERRQ(ierr); 861 ierr = VecGetDM(Grad, &dmGrad);CHKERRQ(ierr); 862 ierr = VecGetArrayRead(Grad, &grad);CHKERRQ(ierr); 863 ierr = PetscFVGetNumComponents(fv, &pdim);CHKERRQ(ierr); 864 ierr = DMGetWorkArray(dm, pdim, MPIU_SCALAR, &fx);CHKERRQ(ierr); 865 } 866 ierr = VecGetArray(locX, &x);CHKERRQ(ierr); 867 for (i = 0; i < numids; ++i) { 868 IS faceIS; 869 const PetscInt *faces; 870 PetscInt numFaces, f; 871 872 ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr); 873 if (!faceIS) continue; /* No points with that id on this process */ 874 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 875 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 876 for (f = 0; f < numFaces; ++f) { 877 const PetscInt face = faces[f], *cells; 878 PetscFVFaceGeom *fg; 879 880 if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */ 881 ierr = PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);CHKERRQ(ierr); 882 if (loc >= 0) continue; 883 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 884 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 885 if (Grad) { 886 PetscFVCellGeom *cg; 887 PetscScalar *cx, *cgrad; 888 PetscScalar *xG; 889 PetscReal dx[3]; 890 PetscInt d; 891 892 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);CHKERRQ(ierr); 893 ierr = DMPlexPointLocalRead(dm, cells[0], x, &cx);CHKERRQ(ierr); 894 ierr = DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);CHKERRQ(ierr); 895 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 896 DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx); 897 for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx); 898 ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx); 899 if (ierru) { 900 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 901 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 902 goto cleanup; 903 } 904 } else { 905 PetscScalar *xI; 906 PetscScalar *xG; 907 908 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr); 909 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 910 ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx); 911 if (ierru) { 912 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 913 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 914 goto cleanup; 915 } 916 } 917 } 918 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 919 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 920 } 921 cleanup: 922 ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr); 923 if (Grad) { 924 ierr = DMRestoreWorkArray(dm, pdim, MPIU_SCALAR, &fx);CHKERRQ(ierr); 925 ierr = VecRestoreArrayRead(Grad, &grad);CHKERRQ(ierr); 926 } 927 if (cellGeometry) {ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr);} 928 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 929 CHKERRQ(ierru); 930 PetscFunctionReturn(0); 931 } 932 933 PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM) 934 { 935 PetscDS prob; 936 PetscInt numBd, b; 937 PetscErrorCode ierr; 938 939 PetscFunctionBegin; 940 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 941 ierr = PetscDSGetNumBoundary(prob, &numBd);CHKERRQ(ierr); 942 for (b = 0; b < numBd; ++b) { 943 DMBoundaryConditionType type; 944 const char *name, *labelname; 945 DMLabel label; 946 PetscInt field, Nc; 947 const PetscInt *comps; 948 PetscObject obj; 949 PetscClassId id; 950 void (*func)(void); 951 PetscInt numids; 952 const PetscInt *ids; 953 void *ctx; 954 955 ierr = DMGetBoundary(dm, b, &type, &name, &labelname, &field, &Nc, &comps, &func, &numids, &ids, &ctx);CHKERRQ(ierr); 956 if (insertEssential != (type & DM_BC_ESSENTIAL)) continue; 957 ierr = DMGetLabel(dm, labelname, &label);CHKERRQ(ierr); 958 if (!label) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "Label %s for boundary condition %s does not exist in the DM", labelname, name); 959 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 960 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 961 if (id == PETSCFE_CLASSID) { 962 switch (type) { 963 /* for FEM, there is no insertion to be done for non-essential boundary conditions */ 964 case DM_BC_ESSENTIAL: 965 ierr = DMPlexLabelAddCells(dm,label);CHKERRQ(ierr); 966 ierr = DMPlexInsertBoundaryValuesEssential(dm, time, field, Nc, comps, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);CHKERRQ(ierr); 967 ierr = DMPlexLabelClearCells(dm,label);CHKERRQ(ierr); 968 break; 969 case DM_BC_ESSENTIAL_FIELD: 970 ierr = DMPlexLabelAddCells(dm,label);CHKERRQ(ierr); 971 ierr = DMPlexInsertBoundaryValuesEssentialField(dm, time, locX, field, Nc, comps, label, numids, ids, 972 (void (*)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 973 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 974 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) func, ctx, locX);CHKERRQ(ierr); 975 ierr = DMPlexLabelClearCells(dm,label);CHKERRQ(ierr); 976 break; 977 default: break; 978 } 979 } else if (id == PETSCFV_CLASSID) { 980 if (!faceGeomFVM) continue; 981 ierr = DMPlexInsertBoundaryValuesRiemann(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, field, Nc, comps, label, numids, ids, 982 (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);CHKERRQ(ierr); 983 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 984 } 985 PetscFunctionReturn(0); 986 } 987 988 /*@ 989 DMPlexInsertBoundaryValues - Puts coefficients which represent boundary values into the local solution vector 990 991 Input Parameters: 992 + dm - The DM 993 . insertEssential - Should I insert essential (e.g. Dirichlet) or inessential (e.g. Neumann) boundary conditions 994 . time - The time 995 . faceGeomFVM - Face geometry data for FV discretizations 996 . cellGeomFVM - Cell geometry data for FV discretizations 997 - gradFVM - Gradient reconstruction data for FV discretizations 998 999 Output Parameters: 1000 . locX - Solution updated with boundary values 1001 1002 Level: developer 1003 1004 .seealso: DMProjectFunctionLabelLocal() 1005 @*/ 1006 PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM) 1007 { 1008 PetscErrorCode ierr; 1009 1010 PetscFunctionBegin; 1011 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1012 PetscValidHeaderSpecific(locX, VEC_CLASSID, 2); 1013 if (faceGeomFVM) {PetscValidHeaderSpecific(faceGeomFVM, VEC_CLASSID, 4);} 1014 if (cellGeomFVM) {PetscValidHeaderSpecific(cellGeomFVM, VEC_CLASSID, 5);} 1015 if (gradFVM) {PetscValidHeaderSpecific(gradFVM, VEC_CLASSID, 6);} 1016 ierr = PetscTryMethod(dm,"DMPlexInsertBoundaryValues_C",(DM,PetscBool,Vec,PetscReal,Vec,Vec,Vec),(dm,insertEssential,locX,time,faceGeomFVM,cellGeomFVM,gradFVM));CHKERRQ(ierr); 1017 PetscFunctionReturn(0); 1018 } 1019 1020 PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 1021 { 1022 Vec localX; 1023 PetscErrorCode ierr; 1024 1025 PetscFunctionBegin; 1026 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1027 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, time, NULL, NULL, NULL);CHKERRQ(ierr); 1028 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1029 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1030 ierr = DMPlexComputeL2DiffLocal(dm, time, funcs, ctxs, localX, diff);CHKERRQ(ierr); 1031 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1032 PetscFunctionReturn(0); 1033 } 1034 1035 /*@C 1036 DMComputeL2DiffLocal - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h. 1037 1038 Collective on DM 1039 1040 Input Parameters: 1041 + dm - The DM 1042 . time - The time 1043 . funcs - The functions to evaluate for each field component 1044 . ctxs - Optional array of contexts to pass to each function, or NULL. 1045 - localX - The coefficient vector u_h, a local vector 1046 1047 Output Parameter: 1048 . diff - The diff ||u - u_h||_2 1049 1050 Level: developer 1051 1052 .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff() 1053 @*/ 1054 PetscErrorCode DMPlexComputeL2DiffLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec localX, PetscReal *diff) 1055 { 1056 const PetscInt debug = ((DM_Plex*)dm->data)->printL2; 1057 DM tdm; 1058 Vec tv; 1059 PetscSection section; 1060 PetscQuadrature quad; 1061 PetscFEGeom fegeom; 1062 PetscScalar *funcVal, *interpolant; 1063 PetscReal *coords, *gcoords; 1064 PetscReal localDiff = 0.0; 1065 const PetscReal *quadWeights; 1066 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cellHeight, cStart, cEnd, cEndInterior, c, field, fieldOffset; 1067 PetscBool transform; 1068 PetscErrorCode ierr; 1069 1070 PetscFunctionBegin; 1071 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1072 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 1073 fegeom.dimEmbed = coordDim; 1074 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 1075 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1076 ierr = DMGetBasisTransformDM_Internal(dm, &tdm);CHKERRQ(ierr); 1077 ierr = DMGetBasisTransformVec_Internal(dm, &tv);CHKERRQ(ierr); 1078 ierr = DMHasBasisTransform(dm, &transform);CHKERRQ(ierr); 1079 for (field = 0; field < numFields; ++field) { 1080 PetscObject obj; 1081 PetscClassId id; 1082 PetscInt Nc; 1083 1084 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 1085 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1086 if (id == PETSCFE_CLASSID) { 1087 PetscFE fe = (PetscFE) obj; 1088 1089 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1090 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1091 } else if (id == PETSCFV_CLASSID) { 1092 PetscFV fv = (PetscFV) obj; 1093 1094 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1095 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1096 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1097 numComponents += Nc; 1098 } 1099 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);CHKERRQ(ierr); 1100 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 1101 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);CHKERRQ(ierr); 1102 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr); 1103 ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr); 1104 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1105 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1106 for (c = cStart; c < cEnd; ++c) { 1107 PetscScalar *x = NULL; 1108 PetscReal elemDiff = 0.0; 1109 PetscInt qc = 0; 1110 1111 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);CHKERRQ(ierr); 1112 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1113 1114 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1115 PetscObject obj; 1116 PetscClassId id; 1117 void * const ctx = ctxs ? ctxs[field] : NULL; 1118 PetscInt Nb, Nc, q, fc; 1119 1120 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 1121 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1122 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1123 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1124 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1125 if (debug) { 1126 char title[1024]; 1127 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 1128 ierr = DMPrintCellVector(c, title, Nb, &x[fieldOffset]);CHKERRQ(ierr); 1129 } 1130 for (q = 0; q < Nq; ++q) { 1131 PetscFEGeom qgeom; 1132 1133 qgeom.dimEmbed = fegeom.dimEmbed; 1134 qgeom.J = &fegeom.J[q*coordDim*coordDim]; 1135 qgeom.invJ = &fegeom.invJ[q*coordDim*coordDim]; 1136 qgeom.detJ = &fegeom.detJ[q]; 1137 if (fegeom.detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, point %D", fegeom.detJ[q], c, q); 1138 if (transform) { 1139 gcoords = &coords[coordDim*Nq]; 1140 ierr = DMPlexBasisTransformApplyReal_Internal(dm, &coords[coordDim*q], PETSC_TRUE, coordDim, &coords[coordDim*q], gcoords, dm->transformCtx);CHKERRQ(ierr); 1141 } else { 1142 gcoords = &coords[coordDim*q]; 1143 } 1144 ierr = (*funcs[field])(coordDim, time, gcoords, Nc, funcVal, ctx); 1145 if (ierr) { 1146 PetscErrorCode ierr2; 1147 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 1148 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 1149 ierr2 = PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2); 1150 CHKERRQ(ierr); 1151 } 1152 if (transform) {ierr = DMPlexBasisTransformApply_Internal(dm, &coords[coordDim*q], PETSC_FALSE, Nc, funcVal, funcVal, dm->transformCtx);CHKERRQ(ierr);} 1153 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], &qgeom, q, interpolant);CHKERRQ(ierr);} 1154 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1155 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1156 for (fc = 0; fc < Nc; ++fc) { 1157 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 1158 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %D field %D,%D diff %g\n", c, field, fc, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q]);CHKERRQ(ierr);} 1159 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q]; 1160 } 1161 } 1162 fieldOffset += Nb; 1163 qc += Nc; 1164 } 1165 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1166 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 1167 localDiff += elemDiff; 1168 } 1169 ierr = PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr); 1170 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1171 *diff = PetscSqrtReal(*diff); 1172 PetscFunctionReturn(0); 1173 } 1174 1175 PetscErrorCode DMComputeL2GradientDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff) 1176 { 1177 const PetscInt debug = ((DM_Plex*)dm->data)->printL2; 1178 DM tdm; 1179 PetscSection section; 1180 PetscQuadrature quad; 1181 Vec localX, tv; 1182 PetscScalar *funcVal, *interpolant; 1183 const PetscReal *quadWeights; 1184 PetscFEGeom fegeom; 1185 PetscReal *coords, *gcoords; 1186 PetscReal localDiff = 0.0; 1187 PetscInt dim, coordDim, qNc = 0, Nq = 0, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset; 1188 PetscBool transform; 1189 PetscErrorCode ierr; 1190 1191 PetscFunctionBegin; 1192 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1193 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 1194 fegeom.dimEmbed = coordDim; 1195 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 1196 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1197 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1198 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1199 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1200 ierr = DMGetBasisTransformDM_Internal(dm, &tdm);CHKERRQ(ierr); 1201 ierr = DMGetBasisTransformVec_Internal(dm, &tv);CHKERRQ(ierr); 1202 ierr = DMHasBasisTransform(dm, &transform);CHKERRQ(ierr); 1203 for (field = 0; field < numFields; ++field) { 1204 PetscFE fe; 1205 PetscInt Nc; 1206 1207 ierr = DMGetField(dm, field, NULL, (PetscObject *) &fe);CHKERRQ(ierr); 1208 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1209 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1210 numComponents += Nc; 1211 } 1212 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);CHKERRQ(ierr); 1213 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 1214 /* ierr = DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 1215 ierr = PetscMalloc6(numComponents,&funcVal,coordDim*Nq,&coords,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ,numComponents*coordDim,&interpolant,Nq,&fegeom.detJ);CHKERRQ(ierr); 1216 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1217 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1218 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1219 for (c = cStart; c < cEnd; ++c) { 1220 PetscScalar *x = NULL; 1221 PetscReal elemDiff = 0.0; 1222 PetscInt qc = 0; 1223 1224 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);CHKERRQ(ierr); 1225 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1226 1227 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1228 PetscFE fe; 1229 void * const ctx = ctxs ? ctxs[field] : NULL; 1230 PetscInt Nb, Nc, q, fc; 1231 1232 ierr = DMGetField(dm, field, NULL, (PetscObject *) &fe);CHKERRQ(ierr); 1233 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1234 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1235 if (debug) { 1236 char title[1024]; 1237 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 1238 ierr = DMPrintCellVector(c, title, Nb, &x[fieldOffset]);CHKERRQ(ierr); 1239 } 1240 for (q = 0; q < Nq; ++q) { 1241 PetscFEGeom qgeom; 1242 1243 qgeom.dimEmbed = fegeom.dimEmbed; 1244 qgeom.J = &fegeom.J[q*coordDim*coordDim]; 1245 qgeom.invJ = &fegeom.invJ[q*coordDim*coordDim]; 1246 qgeom.detJ = &fegeom.detJ[q]; 1247 if (fegeom.detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", fegeom.detJ[q], c, q); 1248 if (transform) { 1249 gcoords = &coords[coordDim*Nq]; 1250 ierr = DMPlexBasisTransformApplyReal_Internal(dm, &coords[coordDim*q], PETSC_TRUE, coordDim, &coords[coordDim*q], gcoords, dm->transformCtx);CHKERRQ(ierr); 1251 } else { 1252 gcoords = &coords[coordDim*q]; 1253 } 1254 ierr = (*funcs[field])(coordDim, time, gcoords, n, Nc, funcVal, ctx); 1255 if (ierr) { 1256 PetscErrorCode ierr2; 1257 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 1258 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 1259 ierr2 = PetscFree6(funcVal,coords,fegeom.J,fegeom.invJ,interpolant,fegeom.detJ);CHKERRQ(ierr2); 1260 CHKERRQ(ierr); 1261 } 1262 if (transform) {ierr = DMPlexBasisTransformApply_Internal(dm, &coords[coordDim*q], PETSC_FALSE, Nc, funcVal, funcVal, dm->transformCtx);CHKERRQ(ierr);} 1263 ierr = PetscFEInterpolateGradient_Static(fe, &x[fieldOffset], &qgeom, q, interpolant);CHKERRQ(ierr); 1264 /* Overwrite with the dot product if the normal is given */ 1265 if (n) { 1266 for (fc = 0; fc < Nc; ++fc) { 1267 PetscScalar sum = 0.0; 1268 PetscInt d; 1269 for (d = 0; d < dim; ++d) sum += interpolant[fc*dim+d]*n[d]; 1270 interpolant[fc] = sum; 1271 } 1272 } 1273 for (fc = 0; fc < Nc; ++fc) { 1274 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 1275 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %D fieldDer %D,%D diff %g\n", c, field, fc, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q]);CHKERRQ(ierr);} 1276 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q]; 1277 } 1278 } 1279 fieldOffset += Nb; 1280 qc += Nc; 1281 } 1282 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1283 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 1284 localDiff += elemDiff; 1285 } 1286 ierr = PetscFree6(funcVal,coords,fegeom.J,fegeom.invJ,interpolant,fegeom.detJ);CHKERRQ(ierr); 1287 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1288 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1289 *diff = PetscSqrtReal(*diff); 1290 PetscFunctionReturn(0); 1291 } 1292 1293 PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 1294 { 1295 const PetscInt debug = ((DM_Plex*)dm->data)->printL2; 1296 DM tdm; 1297 PetscSection section; 1298 PetscQuadrature quad; 1299 Vec localX, tv; 1300 PetscFEGeom fegeom; 1301 PetscScalar *funcVal, *interpolant; 1302 PetscReal *coords, *gcoords; 1303 PetscReal *localDiff; 1304 const PetscReal *quadPoints, *quadWeights; 1305 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset; 1306 PetscBool transform; 1307 PetscErrorCode ierr; 1308 1309 PetscFunctionBegin; 1310 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1311 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 1312 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 1313 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1314 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1315 ierr = VecSet(localX, 0.0);CHKERRQ(ierr); 1316 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1317 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1318 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 1319 ierr = DMGetBasisTransformDM_Internal(dm, &tdm);CHKERRQ(ierr); 1320 ierr = DMGetBasisTransformVec_Internal(dm, &tv);CHKERRQ(ierr); 1321 ierr = DMHasBasisTransform(dm, &transform);CHKERRQ(ierr); 1322 for (field = 0; field < numFields; ++field) { 1323 PetscObject obj; 1324 PetscClassId id; 1325 PetscInt Nc; 1326 1327 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 1328 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1329 if (id == PETSCFE_CLASSID) { 1330 PetscFE fe = (PetscFE) obj; 1331 1332 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1333 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1334 } else if (id == PETSCFV_CLASSID) { 1335 PetscFV fv = (PetscFV) obj; 1336 1337 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1338 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1339 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1340 numComponents += Nc; 1341 } 1342 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 1343 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 1344 ierr = PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,coordDim*(Nq+1),&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);CHKERRQ(ierr); 1345 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1346 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1347 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1348 for (c = cStart; c < cEnd; ++c) { 1349 PetscScalar *x = NULL; 1350 PetscInt qc = 0; 1351 1352 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);CHKERRQ(ierr); 1353 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1354 1355 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1356 PetscObject obj; 1357 PetscClassId id; 1358 void * const ctx = ctxs ? ctxs[field] : NULL; 1359 PetscInt Nb, Nc, q, fc; 1360 1361 PetscReal elemDiff = 0.0; 1362 1363 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 1364 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1365 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1366 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1367 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1368 if (debug) { 1369 char title[1024]; 1370 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 1371 ierr = DMPrintCellVector(c, title, Nb, &x[fieldOffset]);CHKERRQ(ierr); 1372 } 1373 for (q = 0; q < Nq; ++q) { 1374 PetscFEGeom qgeom; 1375 1376 qgeom.dimEmbed = fegeom.dimEmbed; 1377 qgeom.J = &fegeom.J[q*coordDim*coordDim]; 1378 qgeom.invJ = &fegeom.invJ[q*coordDim*coordDim]; 1379 qgeom.detJ = &fegeom.detJ[q]; 1380 if (fegeom.detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature point %D", fegeom.detJ, c, q); 1381 if (transform) { 1382 gcoords = &coords[coordDim*Nq]; 1383 ierr = DMPlexBasisTransformApplyReal_Internal(dm, &coords[coordDim*q], PETSC_TRUE, coordDim, &coords[coordDim*q], gcoords, dm->transformCtx);CHKERRQ(ierr); 1384 } else { 1385 gcoords = &coords[coordDim*q]; 1386 } 1387 ierr = (*funcs[field])(coordDim, time, gcoords, Nc, funcVal, ctx); 1388 if (ierr) { 1389 PetscErrorCode ierr2; 1390 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 1391 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 1392 ierr2 = PetscFree7(localDiff,funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2); 1393 CHKERRQ(ierr); 1394 } 1395 if (transform) {ierr = DMPlexBasisTransformApply_Internal(dm, &coords[coordDim*q], PETSC_FALSE, Nc, funcVal, funcVal, dm->transformCtx);CHKERRQ(ierr);} 1396 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], &qgeom, q, interpolant);CHKERRQ(ierr);} 1397 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1398 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1399 for (fc = 0; fc < Nc; ++fc) { 1400 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 1401 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %D field %D,%D point %g %g %g diff %g\n", c, field, fc, coordDim > 0 ? coords[coordDim*q] : 0., coordDim > 1 ? coords[coordDim*q+1] : 0., coordDim > 2 ? coords[coordDim*q+2] : 0., PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q]);CHKERRQ(ierr);} 1402 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q]; 1403 } 1404 } 1405 fieldOffset += Nb; 1406 qc += Nc; 1407 localDiff[field] += elemDiff; 1408 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d field %d cum diff %g\n", c, field, localDiff[field]);CHKERRQ(ierr);} 1409 } 1410 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1411 } 1412 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1413 ierr = MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1414 for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]); 1415 ierr = PetscFree7(localDiff,funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr); 1416 PetscFunctionReturn(0); 1417 } 1418 1419 /*@C 1420 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. 1421 1422 Collective on DM 1423 1424 Input Parameters: 1425 + dm - The DM 1426 . time - The time 1427 . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation 1428 . ctxs - Optional array of contexts to pass to each function, or NULL. 1429 - X - The coefficient vector u_h 1430 1431 Output Parameter: 1432 . D - A Vec which holds the difference ||u - u_h||_2 for each cell 1433 1434 Level: developer 1435 1436 .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff() 1437 @*/ 1438 PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D) 1439 { 1440 PetscSection section; 1441 PetscQuadrature quad; 1442 Vec localX; 1443 PetscFEGeom fegeom; 1444 PetscScalar *funcVal, *interpolant; 1445 PetscReal *coords; 1446 const PetscReal *quadPoints, *quadWeights; 1447 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset; 1448 PetscErrorCode ierr; 1449 1450 PetscFunctionBegin; 1451 ierr = VecSet(D, 0.0);CHKERRQ(ierr); 1452 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1453 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 1454 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 1455 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1456 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1457 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 1458 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1459 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1460 for (field = 0; field < numFields; ++field) { 1461 PetscObject obj; 1462 PetscClassId id; 1463 PetscInt Nc; 1464 1465 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 1466 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1467 if (id == PETSCFE_CLASSID) { 1468 PetscFE fe = (PetscFE) obj; 1469 1470 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1471 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1472 } else if (id == PETSCFV_CLASSID) { 1473 PetscFV fv = (PetscFV) obj; 1474 1475 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1476 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1477 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1478 numComponents += Nc; 1479 } 1480 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 1481 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 1482 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);CHKERRQ(ierr); 1483 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1484 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1485 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1486 for (c = cStart; c < cEnd; ++c) { 1487 PetscScalar *x = NULL; 1488 PetscScalar elemDiff = 0.0; 1489 PetscInt qc = 0; 1490 1491 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);CHKERRQ(ierr); 1492 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1493 1494 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1495 PetscObject obj; 1496 PetscClassId id; 1497 void * const ctx = ctxs ? ctxs[field] : NULL; 1498 PetscInt Nb, Nc, q, fc; 1499 1500 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 1501 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1502 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1503 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1504 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1505 if (funcs[field]) { 1506 for (q = 0; q < Nq; ++q) { 1507 PetscFEGeom qgeom; 1508 1509 qgeom.dimEmbed = fegeom.dimEmbed; 1510 qgeom.J = &fegeom.J[q*coordDim*coordDim]; 1511 qgeom.invJ = &fegeom.invJ[q*coordDim*coordDim]; 1512 qgeom.detJ = &fegeom.detJ[q]; 1513 if (fegeom.detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", (double)fegeom.detJ[q], c, q); 1514 ierr = (*funcs[field])(coordDim, time, &coords[q*coordDim], Nc, funcVal, ctx); 1515 if (ierr) { 1516 PetscErrorCode ierr2; 1517 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 1518 ierr2 = PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2); 1519 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 1520 CHKERRQ(ierr); 1521 } 1522 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], &qgeom, q, interpolant);CHKERRQ(ierr);} 1523 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1524 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1525 for (fc = 0; fc < Nc; ++fc) { 1526 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 1527 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q]; 1528 } 1529 } 1530 } 1531 fieldOffset += Nb; 1532 qc += Nc; 1533 } 1534 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1535 ierr = VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);CHKERRQ(ierr); 1536 } 1537 ierr = PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr); 1538 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1539 ierr = VecSqrtAbs(D);CHKERRQ(ierr); 1540 PetscFunctionReturn(0); 1541 } 1542 1543 /*@C 1544 DMPlexComputeGradientClementInterpolant - This function computes the L2 projection of the cellwise gradient of a function u onto P1, and stores it in a Vec. 1545 1546 Collective on DM 1547 1548 Input Parameters: 1549 + dm - The DM 1550 - LocX - The coefficient vector u_h 1551 1552 Output Parameter: 1553 . locC - A Vec which holds the Clement interpolant of the gradient 1554 1555 Notes: 1556 Add citation to (Clement, 1975) and definition of the interpolant 1557 \nabla u_h(v_i) = \sum_{T_i \in support(v_i)} |T_i| \nabla u_h(T_i) / \sum_{T_i \in support(v_i)} |T_i| where |T_i| is the cell volume 1558 1559 Level: developer 1560 1561 .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff() 1562 @*/ 1563 PetscErrorCode DMPlexComputeGradientClementInterpolant(DM dm, Vec locX, Vec locC) 1564 { 1565 DM_Plex *mesh = (DM_Plex *) dm->data; 1566 PetscInt debug = mesh->printFEM; 1567 DM dmC; 1568 PetscSection section; 1569 PetscQuadrature quad; 1570 PetscScalar *interpolant, *gradsum; 1571 PetscFEGeom fegeom; 1572 PetscReal *coords; 1573 const PetscReal *quadPoints, *quadWeights; 1574 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, vStart, vEnd, v, field, fieldOffset; 1575 PetscErrorCode ierr; 1576 1577 PetscFunctionBegin; 1578 ierr = VecGetDM(locC, &dmC);CHKERRQ(ierr); 1579 ierr = VecSet(locC, 0.0);CHKERRQ(ierr); 1580 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1581 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 1582 fegeom.dimEmbed = coordDim; 1583 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 1584 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1585 for (field = 0; field < numFields; ++field) { 1586 PetscObject obj; 1587 PetscClassId id; 1588 PetscInt Nc; 1589 1590 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 1591 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1592 if (id == PETSCFE_CLASSID) { 1593 PetscFE fe = (PetscFE) obj; 1594 1595 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1596 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1597 } else if (id == PETSCFV_CLASSID) { 1598 PetscFV fv = (PetscFV) obj; 1599 1600 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1601 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1602 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1603 numComponents += Nc; 1604 } 1605 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 1606 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 1607 ierr = PetscMalloc6(coordDim*numComponents*2,&gradsum,coordDim*numComponents,&interpolant,coordDim*Nq,&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);CHKERRQ(ierr); 1608 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1609 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1610 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1611 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1612 for (v = vStart; v < vEnd; ++v) { 1613 PetscScalar volsum = 0.0; 1614 PetscInt *star = NULL; 1615 PetscInt starSize, st, d, fc; 1616 1617 ierr = PetscMemzero(gradsum, coordDim*numComponents * sizeof(PetscScalar));CHKERRQ(ierr); 1618 ierr = DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1619 for (st = 0; st < starSize*2; st += 2) { 1620 const PetscInt cell = star[st]; 1621 PetscScalar *grad = &gradsum[coordDim*numComponents]; 1622 PetscScalar *x = NULL; 1623 PetscReal vol = 0.0; 1624 1625 if ((cell < cStart) || (cell >= cEnd)) continue; 1626 ierr = DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);CHKERRQ(ierr); 1627 ierr = DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr); 1628 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1629 PetscObject obj; 1630 PetscClassId id; 1631 PetscInt Nb, Nc, q, qc = 0; 1632 1633 ierr = PetscMemzero(grad, coordDim*numComponents * sizeof(PetscScalar));CHKERRQ(ierr); 1634 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 1635 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1636 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1637 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1638 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1639 for (q = 0; q < Nq; ++q) { 1640 PetscFEGeom qgeom; 1641 1642 qgeom.dimEmbed = fegeom.dimEmbed; 1643 qgeom.J = &fegeom.J[q*coordDim*coordDim]; 1644 qgeom.invJ = &fegeom.invJ[q*coordDim*coordDim]; 1645 qgeom.detJ = &fegeom.detJ[q]; 1646 if (fegeom.detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", (double)fegeom.detJ[q], cell, q); 1647 if (ierr) { 1648 PetscErrorCode ierr2; 1649 ierr2 = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr2); 1650 ierr2 = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr2); 1651 ierr2 = PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2); 1652 CHKERRQ(ierr); 1653 } 1654 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolateGradient_Static((PetscFE) obj, &x[fieldOffset], &qgeom, q, interpolant);CHKERRQ(ierr);} 1655 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 1656 for (fc = 0; fc < Nc; ++fc) { 1657 const PetscReal wt = quadWeights[q*qNc+qc+fc]; 1658 1659 for (d = 0; d < coordDim; ++d) grad[fc*coordDim+d] += interpolant[fc*dim+d]*wt*fegeom.detJ[q]; 1660 } 1661 vol += quadWeights[q*qNc]*fegeom.detJ[q]; 1662 } 1663 fieldOffset += Nb; 1664 qc += Nc; 1665 } 1666 ierr = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr); 1667 for (fc = 0; fc < numComponents; ++fc) { 1668 for (d = 0; d < coordDim; ++d) { 1669 gradsum[fc*coordDim+d] += grad[fc*coordDim+d]; 1670 } 1671 } 1672 volsum += vol; 1673 if (debug) { 1674 ierr = PetscPrintf(PETSC_COMM_SELF, "Cell %D gradient: [", cell);CHKERRQ(ierr); 1675 for (fc = 0; fc < numComponents; ++fc) { 1676 for (d = 0; d < coordDim; ++d) { 1677 if (fc || d > 0) {ierr = PetscPrintf(PETSC_COMM_SELF, ", ");CHKERRQ(ierr);} 1678 ierr = PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc*coordDim+d]));CHKERRQ(ierr); 1679 } 1680 } 1681 ierr = PetscPrintf(PETSC_COMM_SELF, "]\n");CHKERRQ(ierr); 1682 } 1683 } 1684 for (fc = 0; fc < numComponents; ++fc) { 1685 for (d = 0; d < coordDim; ++d) gradsum[fc*coordDim+d] /= volsum; 1686 } 1687 ierr = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1688 ierr = DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES);CHKERRQ(ierr); 1689 } 1690 ierr = PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr); 1691 PetscFunctionReturn(0); 1692 } 1693 1694 static PetscErrorCode DMPlexComputeIntegral_Internal(DM dm, Vec X, PetscInt cStart, PetscInt cEnd, PetscScalar *cintegral, void *user) 1695 { 1696 DM dmAux = NULL; 1697 PetscDS prob, probAux = NULL; 1698 PetscSection section, sectionAux; 1699 Vec locX, locA; 1700 PetscInt dim, numCells = cEnd - cStart, c, f; 1701 PetscBool useFVM = PETSC_FALSE; 1702 /* DS */ 1703 PetscInt Nf, totDim, *uOff, *uOff_x, numConstants; 1704 PetscInt NfAux, totDimAux, *aOff; 1705 PetscScalar *u, *a; 1706 const PetscScalar *constants; 1707 /* Geometry */ 1708 PetscFEGeom *cgeomFEM; 1709 DM dmGrad; 1710 PetscQuadrature affineQuad = NULL; 1711 Vec cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL; 1712 PetscFVCellGeom *cgeomFVM; 1713 const PetscScalar *lgrad; 1714 PetscInt maxDegree; 1715 DMField coordField; 1716 IS cellIS; 1717 PetscErrorCode ierr; 1718 1719 PetscFunctionBegin; 1720 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1721 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1722 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 1723 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1724 /* Determine which discretizations we have */ 1725 for (f = 0; f < Nf; ++f) { 1726 PetscObject obj; 1727 PetscClassId id; 1728 1729 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1730 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1731 if (id == PETSCFV_CLASSID) useFVM = PETSC_TRUE; 1732 } 1733 /* Get local solution with boundary values */ 1734 ierr = DMGetLocalVector(dm, &locX);CHKERRQ(ierr); 1735 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 1736 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);CHKERRQ(ierr); 1737 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);CHKERRQ(ierr); 1738 /* Read DS information */ 1739 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1740 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 1741 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 1742 ierr = ISCreateStride(PETSC_COMM_SELF,numCells,cStart,1,&cellIS);CHKERRQ(ierr); 1743 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 1744 /* Read Auxiliary DS information */ 1745 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1746 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);CHKERRQ(ierr); 1747 if (dmAux) { 1748 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1749 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 1750 ierr = DMGetSection(dmAux, §ionAux);CHKERRQ(ierr); 1751 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1752 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 1753 } 1754 /* Allocate data arrays */ 1755 ierr = PetscCalloc1(numCells*totDim, &u);CHKERRQ(ierr); 1756 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1757 /* Read out geometry */ 1758 ierr = DMGetCoordinateField(dm,&coordField);CHKERRQ(ierr); 1759 ierr = DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);CHKERRQ(ierr); 1760 if (maxDegree <= 1) { 1761 ierr = DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);CHKERRQ(ierr); 1762 if (affineQuad) { 1763 ierr = DMFieldCreateFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&cgeomFEM);CHKERRQ(ierr); 1764 } 1765 } 1766 if (useFVM) { 1767 PetscFV fv = NULL; 1768 Vec grad; 1769 PetscInt fStart, fEnd; 1770 PetscBool compGrad; 1771 1772 for (f = 0; f < Nf; ++f) { 1773 PetscObject obj; 1774 PetscClassId id; 1775 1776 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1777 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1778 if (id == PETSCFV_CLASSID) {fv = (PetscFV) obj; break;} 1779 } 1780 ierr = PetscFVGetComputeGradients(fv, &compGrad);CHKERRQ(ierr); 1781 ierr = PetscFVSetComputeGradients(fv, PETSC_TRUE);CHKERRQ(ierr); 1782 ierr = DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);CHKERRQ(ierr); 1783 ierr = DMPlexComputeGradientFVM(dm, fv, faceGeometryFVM, cellGeometryFVM, &dmGrad);CHKERRQ(ierr); 1784 ierr = PetscFVSetComputeGradients(fv, compGrad);CHKERRQ(ierr); 1785 ierr = VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1786 /* Reconstruct and limit cell gradients */ 1787 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1788 ierr = DMGetGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1789 ierr = DMPlexReconstructGradients_Internal(dm, fv, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);CHKERRQ(ierr); 1790 /* Communicate gradient values */ 1791 ierr = DMGetLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1792 ierr = DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1793 ierr = DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1794 ierr = DMRestoreGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1795 /* Handle non-essential (e.g. outflow) boundary values */ 1796 ierr = DMPlexInsertBoundaryValues(dm, PETSC_FALSE, locX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);CHKERRQ(ierr); 1797 ierr = VecGetArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1798 } 1799 /* Read out data from inputs */ 1800 for (c = cStart; c < cEnd; ++c) { 1801 PetscScalar *x = NULL; 1802 PetscInt i; 1803 1804 ierr = DMPlexVecGetClosure(dm, section, locX, c, NULL, &x);CHKERRQ(ierr); 1805 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 1806 ierr = DMPlexVecRestoreClosure(dm, section, locX, c, NULL, &x);CHKERRQ(ierr); 1807 if (dmAux) { 1808 ierr = DMPlexVecGetClosure(dmAux, sectionAux, locA, c, NULL, &x);CHKERRQ(ierr); 1809 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 1810 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, locA, c, NULL, &x);CHKERRQ(ierr); 1811 } 1812 } 1813 /* Do integration for each field */ 1814 for (f = 0; f < Nf; ++f) { 1815 PetscObject obj; 1816 PetscClassId id; 1817 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 1818 1819 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1820 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1821 if (id == PETSCFE_CLASSID) { 1822 PetscFE fe = (PetscFE) obj; 1823 PetscQuadrature q; 1824 PetscFEGeom *chunkGeom = NULL; 1825 PetscInt Nq, Nb; 1826 1827 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1828 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1829 ierr = PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 1830 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1831 blockSize = Nb*Nq; 1832 batchSize = numBlocks * blockSize; 1833 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1834 numChunks = numCells / (numBatches*batchSize); 1835 Ne = numChunks*numBatches*batchSize; 1836 Nr = numCells % (numBatches*batchSize); 1837 offset = numCells - Nr; 1838 if (!affineQuad) { 1839 ierr = DMFieldCreateFEGeom(coordField,cellIS,q,PETSC_FALSE,&cgeomFEM);CHKERRQ(ierr); 1840 } 1841 ierr = PetscFEGeomGetChunk(cgeomFEM,0,offset,&chunkGeom);CHKERRQ(ierr); 1842 ierr = PetscFEIntegrate(prob, f, Ne, chunkGeom, u, probAux, a, cintegral);CHKERRQ(ierr); 1843 ierr = PetscFEGeomGetChunk(cgeomFEM,offset,numCells,&chunkGeom);CHKERRQ(ierr); 1844 ierr = PetscFEIntegrate(prob, f, Nr, chunkGeom, &u[offset*totDim], probAux, &a[offset*totDimAux], &cintegral[offset*Nf]);CHKERRQ(ierr); 1845 ierr = PetscFEGeomRestoreChunk(cgeomFEM,offset,numCells,&chunkGeom);CHKERRQ(ierr); 1846 if (!affineQuad) { 1847 ierr = PetscFEGeomDestroy(&cgeomFEM);CHKERRQ(ierr); 1848 } 1849 } else if (id == PETSCFV_CLASSID) { 1850 PetscInt foff; 1851 PetscPointFunc obj_func; 1852 PetscScalar lint; 1853 1854 ierr = PetscDSGetObjective(prob, f, &obj_func);CHKERRQ(ierr); 1855 ierr = PetscDSGetFieldOffset(prob, f, &foff);CHKERRQ(ierr); 1856 if (obj_func) { 1857 for (c = 0; c < numCells; ++c) { 1858 PetscScalar *u_x; 1859 1860 ierr = DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);CHKERRQ(ierr); 1861 obj_func(dim, Nf, NfAux, uOff, uOff_x, &u[totDim*c+foff], NULL, u_x, aOff, NULL, &a[totDimAux*c], NULL, NULL, 0.0, cgeomFVM[c].centroid, numConstants, constants, &lint); 1862 cintegral[c*Nf+f] += PetscRealPart(lint)*cgeomFVM[c].volume; 1863 } 1864 } 1865 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 1866 } 1867 /* Cleanup data arrays */ 1868 if (useFVM) { 1869 ierr = VecRestoreArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1870 ierr = VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1871 ierr = DMRestoreLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1872 ierr = VecDestroy(&faceGeometryFVM);CHKERRQ(ierr); 1873 ierr = VecDestroy(&cellGeometryFVM);CHKERRQ(ierr); 1874 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr); 1875 } 1876 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1877 ierr = PetscFree(u);CHKERRQ(ierr); 1878 /* Cleanup */ 1879 if (affineQuad) { 1880 ierr = PetscFEGeomDestroy(&cgeomFEM);CHKERRQ(ierr); 1881 } 1882 ierr = PetscQuadratureDestroy(&affineQuad);CHKERRQ(ierr); 1883 ierr = ISDestroy(&cellIS);CHKERRQ(ierr); 1884 ierr = DMRestoreLocalVector(dm, &locX);CHKERRQ(ierr); 1885 PetscFunctionReturn(0); 1886 } 1887 1888 /*@ 1889 DMPlexComputeIntegralFEM - Form the integral over the domain from the global input X using pointwise functions specified by the user 1890 1891 Input Parameters: 1892 + dm - The mesh 1893 . X - Global input vector 1894 - user - The user context 1895 1896 Output Parameter: 1897 . integral - Integral for each field 1898 1899 Level: developer 1900 1901 .seealso: DMPlexComputeResidualFEM() 1902 @*/ 1903 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscScalar *integral, void *user) 1904 { 1905 DM_Plex *mesh = (DM_Plex *) dm->data; 1906 PetscScalar *cintegral, *lintegral; 1907 PetscInt Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell; 1908 PetscErrorCode ierr; 1909 1910 PetscFunctionBegin; 1911 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1912 PetscValidHeaderSpecific(X, VEC_CLASSID, 2); 1913 PetscValidPointer(integral, 3); 1914 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1915 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 1916 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr); 1917 ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr); 1918 ierr = DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);CHKERRQ(ierr); 1919 cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight]; 1920 /* TODO Introduce a loop over large chunks (right now this is a single chunk) */ 1921 ierr = PetscCalloc2(Nf, &lintegral, (cEnd-cStart)*Nf, &cintegral);CHKERRQ(ierr); 1922 ierr = DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);CHKERRQ(ierr); 1923 /* Sum up values */ 1924 for (cell = cStart; cell < cEnd; ++cell) { 1925 const PetscInt c = cell - cStart; 1926 1927 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);CHKERRQ(ierr);} 1928 for (f = 0; f < Nf; ++f) lintegral[f] += cintegral[c*Nf+f]; 1929 } 1930 ierr = MPIU_Allreduce(lintegral, integral, Nf, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr); 1931 if (mesh->printFEM) { 1932 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Integral:");CHKERRQ(ierr); 1933 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", (double) PetscRealPart(integral[f]));CHKERRQ(ierr);} 1934 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 1935 } 1936 ierr = PetscFree2(lintegral, cintegral);CHKERRQ(ierr); 1937 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1938 PetscFunctionReturn(0); 1939 } 1940 1941 /*@ 1942 DMPlexComputeCellwiseIntegralFEM - Form the vector of cellwise integrals F from the global input X using pointwise functions specified by the user 1943 1944 Input Parameters: 1945 + dm - The mesh 1946 . X - Global input vector 1947 - user - The user context 1948 1949 Output Parameter: 1950 . integral - Cellwise integrals for each field 1951 1952 Level: developer 1953 1954 .seealso: DMPlexComputeResidualFEM() 1955 @*/ 1956 PetscErrorCode DMPlexComputeCellwiseIntegralFEM(DM dm, Vec X, Vec F, void *user) 1957 { 1958 DM_Plex *mesh = (DM_Plex *) dm->data; 1959 DM dmF; 1960 PetscSection sectionF; 1961 PetscScalar *cintegral, *af; 1962 PetscInt Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell; 1963 PetscErrorCode ierr; 1964 1965 PetscFunctionBegin; 1966 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1967 PetscValidHeaderSpecific(X, VEC_CLASSID, 2); 1968 PetscValidHeaderSpecific(F, VEC_CLASSID, 3); 1969 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1970 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 1971 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr); 1972 ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr); 1973 ierr = DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);CHKERRQ(ierr); 1974 cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight]; 1975 /* TODO Introduce a loop over large chunks (right now this is a single chunk) */ 1976 ierr = PetscCalloc1((cEnd-cStart)*Nf, &cintegral);CHKERRQ(ierr); 1977 ierr = DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);CHKERRQ(ierr); 1978 /* Put values in F*/ 1979 ierr = VecGetDM(F, &dmF);CHKERRQ(ierr); 1980 ierr = DMGetSection(dmF, §ionF);CHKERRQ(ierr); 1981 ierr = VecGetArray(F, &af);CHKERRQ(ierr); 1982 for (cell = cStart; cell < cEnd; ++cell) { 1983 const PetscInt c = cell - cStart; 1984 PetscInt dof, off; 1985 1986 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);CHKERRQ(ierr);} 1987 ierr = PetscSectionGetDof(sectionF, cell, &dof);CHKERRQ(ierr); 1988 ierr = PetscSectionGetOffset(sectionF, cell, &off);CHKERRQ(ierr); 1989 if (dof != Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cell dofs %D != %D", dof, Nf); 1990 for (f = 0; f < Nf; ++f) af[off+f] = cintegral[c*Nf+f]; 1991 } 1992 ierr = VecRestoreArray(F, &af);CHKERRQ(ierr); 1993 ierr = PetscFree(cintegral);CHKERRQ(ierr); 1994 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1995 PetscFunctionReturn(0); 1996 } 1997 1998 static PetscErrorCode DMPlexComputeBdIntegral_Internal(DM dm, Vec locX, IS pointIS, 1999 void (*func)(PetscInt, PetscInt, PetscInt, 2000 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 2001 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 2002 PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 2003 PetscScalar *fintegral, void *user) 2004 { 2005 DM plex = NULL, plexA = NULL; 2006 PetscDS prob, probAux = NULL; 2007 PetscSection section, sectionAux = NULL; 2008 Vec locA = NULL; 2009 DMField coordField; 2010 PetscInt Nf, totDim, *uOff, *uOff_x; 2011 PetscInt NfAux = 0, totDimAux = 0, *aOff = NULL; 2012 PetscScalar *u, *a = NULL; 2013 const PetscScalar *constants; 2014 PetscInt numConstants, f; 2015 PetscErrorCode ierr; 2016 2017 PetscFunctionBegin; 2018 ierr = DMGetCoordinateField(dm, &coordField);CHKERRQ(ierr); 2019 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 2020 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 2021 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 2022 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 2023 /* Determine which discretizations we have */ 2024 for (f = 0; f < Nf; ++f) { 2025 PetscObject obj; 2026 PetscClassId id; 2027 2028 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 2029 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2030 if (id == PETSCFV_CLASSID) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Not supported for FVM (field %D)", f); 2031 } 2032 /* Read DS information */ 2033 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 2034 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 2035 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 2036 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 2037 /* Read Auxiliary DS information */ 2038 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);CHKERRQ(ierr); 2039 if (locA) { 2040 DM dmAux; 2041 2042 ierr = VecGetDM(locA, &dmAux);CHKERRQ(ierr); 2043 ierr = DMConvert(dmAux, DMPLEX, &plexA);CHKERRQ(ierr); 2044 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 2045 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 2046 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 2047 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 2048 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 2049 } 2050 /* Integrate over points */ 2051 { 2052 PetscFEGeom *fgeom, *chunkGeom = NULL; 2053 PetscInt maxDegree; 2054 PetscQuadrature qGeom = NULL; 2055 const PetscInt *points; 2056 PetscInt numFaces, face, Nq, field; 2057 PetscInt numChunks, chunkSize, chunk, Nr, offset; 2058 2059 ierr = ISGetLocalSize(pointIS, &numFaces);CHKERRQ(ierr); 2060 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 2061 ierr = PetscCalloc2(numFaces*totDim, &u, locA ? numFaces*totDimAux : 0, &a);CHKERRQ(ierr); 2062 ierr = DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree);CHKERRQ(ierr); 2063 for (field = 0; field < Nf; ++field) { 2064 PetscFE fe; 2065 2066 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 2067 if (maxDegree <= 1) {ierr = DMFieldCreateDefaultQuadrature(coordField, pointIS, &qGeom);CHKERRQ(ierr);} 2068 if (!qGeom) { 2069 ierr = PetscFEGetFaceQuadrature(fe, &qGeom);CHKERRQ(ierr); 2070 ierr = PetscObjectReference((PetscObject) qGeom);CHKERRQ(ierr); 2071 } 2072 ierr = PetscQuadratureGetData(qGeom, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 2073 ierr = DMPlexGetFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom);CHKERRQ(ierr); 2074 for (face = 0; face < numFaces; ++face) { 2075 const PetscInt point = points[face], *support, *cone; 2076 PetscScalar *x = NULL; 2077 PetscInt i, coneSize, faceLoc; 2078 2079 ierr = DMPlexGetSupport(dm, point, &support);CHKERRQ(ierr); 2080 ierr = DMPlexGetConeSize(dm, support[0], &coneSize);CHKERRQ(ierr); 2081 ierr = DMPlexGetCone(dm, support[0], &cone);CHKERRQ(ierr); 2082 for (faceLoc = 0; faceLoc < coneSize; ++faceLoc) if (cone[faceLoc] == point) break; 2083 if (faceLoc == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %D in cone of support[0] %D", face, support[0]); 2084 fgeom->face[face][0] = faceLoc; 2085 ierr = DMPlexVecGetClosure(plex, section, locX, support[0], NULL, &x);CHKERRQ(ierr); 2086 for (i = 0; i < totDim; ++i) u[face*totDim+i] = x[i]; 2087 ierr = DMPlexVecRestoreClosure(plex, section, locX, support[0], NULL, &x);CHKERRQ(ierr); 2088 if (locA) { 2089 PetscInt subp; 2090 ierr = DMPlexGetSubpoint(plexA, support[0], &subp);CHKERRQ(ierr); 2091 ierr = DMPlexVecGetClosure(plexA, sectionAux, locA, subp, NULL, &x);CHKERRQ(ierr); 2092 for (i = 0; i < totDimAux; ++i) a[f*totDimAux+i] = x[i]; 2093 ierr = DMPlexVecRestoreClosure(plexA, sectionAux, locA, subp, NULL, &x);CHKERRQ(ierr); 2094 } 2095 } 2096 /* Get blocking */ 2097 { 2098 PetscQuadrature q; 2099 PetscInt numBatches, batchSize, numBlocks, blockSize; 2100 PetscInt Nq, Nb; 2101 2102 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 2103 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 2104 ierr = PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 2105 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 2106 blockSize = Nb*Nq; 2107 batchSize = numBlocks * blockSize; 2108 chunkSize = numBatches*batchSize; 2109 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 2110 numChunks = numFaces / chunkSize; 2111 Nr = numFaces % chunkSize; 2112 offset = numFaces - Nr; 2113 } 2114 /* Do integration for each field */ 2115 for (chunk = 0; chunk < numChunks; ++chunk) { 2116 ierr = PetscFEGeomGetChunk(fgeom, chunk*chunkSize, (chunk+1)*chunkSize, &chunkGeom);CHKERRQ(ierr); 2117 ierr = PetscFEIntegrateBd(prob, field, func, chunkSize, chunkGeom, u, probAux, a, fintegral);CHKERRQ(ierr); 2118 ierr = PetscFEGeomRestoreChunk(fgeom, 0, offset, &chunkGeom);CHKERRQ(ierr); 2119 } 2120 ierr = PetscFEGeomGetChunk(fgeom, offset, numFaces, &chunkGeom);CHKERRQ(ierr); 2121 ierr = PetscFEIntegrateBd(prob, field, func, Nr, chunkGeom, &u[offset*totDim], probAux, a ? &a[offset*totDimAux] : NULL, &fintegral[offset*Nf]);CHKERRQ(ierr); 2122 ierr = PetscFEGeomRestoreChunk(fgeom, offset, numFaces, &chunkGeom);CHKERRQ(ierr); 2123 /* Cleanup data arrays */ 2124 ierr = DMPlexRestoreFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom);CHKERRQ(ierr); 2125 ierr = PetscQuadratureDestroy(&qGeom);CHKERRQ(ierr); 2126 ierr = PetscFree2(u, a);CHKERRQ(ierr); 2127 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 2128 } 2129 } 2130 if (plex) {ierr = DMDestroy(&plex);CHKERRQ(ierr);} 2131 if (plexA) {ierr = DMDestroy(&plexA);CHKERRQ(ierr);} 2132 PetscFunctionReturn(0); 2133 } 2134 2135 /*@ 2136 DMPlexComputeBdIntegral - Form the integral over the specified boundary from the global input X using pointwise functions specified by the user 2137 2138 Input Parameters: 2139 + dm - The mesh 2140 . X - Global input vector 2141 . label - The boundary DMLabel 2142 . numVals - The number of label values to use, or PETSC_DETERMINE for all values 2143 . vals - The label values to use, or PETSC_NULL for all values 2144 . func = The function to integrate along the boundary 2145 - user - The user context 2146 2147 Output Parameter: 2148 . integral - Integral for each field 2149 2150 Level: developer 2151 2152 .seealso: DMPlexComputeIntegralFEM(), DMPlexComputeBdResidualFEM() 2153 @*/ 2154 PetscErrorCode DMPlexComputeBdIntegral(DM dm, Vec X, DMLabel label, PetscInt numVals, const PetscInt vals[], 2155 void (*func)(PetscInt, PetscInt, PetscInt, 2156 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 2157 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 2158 PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 2159 PetscScalar *integral, void *user) 2160 { 2161 Vec locX; 2162 PetscSection section; 2163 DMLabel depthLabel; 2164 IS facetIS; 2165 PetscInt dim, Nf, f, v; 2166 PetscErrorCode ierr; 2167 2168 PetscFunctionBegin; 2169 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2170 PetscValidHeaderSpecific(X, VEC_CLASSID, 2); 2171 PetscValidPointer(label, 3); 2172 if (vals) PetscValidPointer(vals, 5); 2173 PetscValidPointer(integral, 6); 2174 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 2175 ierr = DMPlexGetDepthLabel(dm, &depthLabel);CHKERRQ(ierr); 2176 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2177 ierr = DMLabelGetStratumIS(depthLabel, dim-1, &facetIS);CHKERRQ(ierr); 2178 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 2179 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 2180 /* Get local solution with boundary values */ 2181 ierr = DMGetLocalVector(dm, &locX);CHKERRQ(ierr); 2182 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 2183 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);CHKERRQ(ierr); 2184 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);CHKERRQ(ierr); 2185 /* Loop over label values */ 2186 ierr = PetscMemzero(integral, Nf * sizeof(PetscScalar));CHKERRQ(ierr); 2187 for (v = 0; v < numVals; ++v) { 2188 IS pointIS; 2189 PetscInt numFaces, face; 2190 PetscScalar *fintegral; 2191 2192 ierr = DMLabelGetStratumIS(label, vals[v], &pointIS);CHKERRQ(ierr); 2193 if (!pointIS) continue; /* No points with that id on this process */ 2194 { 2195 IS isectIS; 2196 2197 /* TODO: Special cases of ISIntersect where it is quick to check a priori if one is a superset of the other */ 2198 ierr = ISIntersect_Caching_Internal(facetIS, pointIS, &isectIS);CHKERRQ(ierr); 2199 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 2200 pointIS = isectIS; 2201 } 2202 ierr = ISGetLocalSize(pointIS, &numFaces);CHKERRQ(ierr); 2203 ierr = PetscCalloc1(numFaces*Nf, &fintegral);CHKERRQ(ierr); 2204 ierr = DMPlexComputeBdIntegral_Internal(dm, locX, pointIS, func, fintegral, user);CHKERRQ(ierr); 2205 /* Sum point contributions into integral */ 2206 for (f = 0; f < Nf; ++f) for (face = 0; face < numFaces; ++face) integral[f] += fintegral[face*Nf+f]; 2207 ierr = PetscFree(fintegral);CHKERRQ(ierr); 2208 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 2209 } 2210 ierr = DMRestoreLocalVector(dm, &locX);CHKERRQ(ierr); 2211 ierr = ISDestroy(&facetIS);CHKERRQ(ierr); 2212 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 2213 PetscFunctionReturn(0); 2214 } 2215 2216 /*@ 2217 DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 2218 2219 Input Parameters: 2220 + dmf - The fine mesh 2221 . dmc - The coarse mesh 2222 - user - The user context 2223 2224 Output Parameter: 2225 . In - The interpolation matrix 2226 2227 Level: developer 2228 2229 .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 2230 @*/ 2231 PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user) 2232 { 2233 DM_Plex *mesh = (DM_Plex *) dmc->data; 2234 const char *name = "Interpolator"; 2235 PetscDS prob; 2236 PetscFE *feRef; 2237 PetscFV *fvRef; 2238 PetscSection fsection, fglobalSection; 2239 PetscSection csection, cglobalSection; 2240 PetscScalar *elemMat; 2241 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c; 2242 PetscInt cTotDim, rTotDim = 0; 2243 PetscErrorCode ierr; 2244 2245 PetscFunctionBegin; 2246 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2247 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 2248 ierr = DMGetSection(dmf, &fsection);CHKERRQ(ierr); 2249 ierr = DMGetGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 2250 ierr = DMGetSection(dmc, &csection);CHKERRQ(ierr); 2251 ierr = DMGetGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 2252 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 2253 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 2254 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 2255 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 2256 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr); 2257 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 2258 for (f = 0; f < Nf; ++f) { 2259 PetscObject obj; 2260 PetscClassId id; 2261 PetscInt rNb = 0, Nc = 0; 2262 2263 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 2264 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2265 if (id == PETSCFE_CLASSID) { 2266 PetscFE fe = (PetscFE) obj; 2267 2268 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 2269 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 2270 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 2271 } else if (id == PETSCFV_CLASSID) { 2272 PetscFV fv = (PetscFV) obj; 2273 PetscDualSpace Q; 2274 2275 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 2276 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 2277 ierr = PetscDualSpaceGetDimension(Q, &rNb);CHKERRQ(ierr); 2278 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 2279 } 2280 rTotDim += rNb; 2281 } 2282 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 2283 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 2284 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 2285 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 2286 PetscDualSpace Qref; 2287 PetscQuadrature f; 2288 const PetscReal *qpoints, *qweights; 2289 PetscReal *points; 2290 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 2291 2292 /* Compose points from all dual basis functionals */ 2293 if (feRef[fieldI]) { 2294 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 2295 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 2296 } else { 2297 ierr = PetscFVGetDualSpace(fvRef[fieldI], &Qref);CHKERRQ(ierr); 2298 ierr = PetscFVGetNumComponents(fvRef[fieldI], &Nc);CHKERRQ(ierr); 2299 } 2300 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 2301 for (i = 0; i < fpdim; ++i) { 2302 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 2303 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 2304 npoints += Np; 2305 } 2306 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 2307 for (i = 0, k = 0; i < fpdim; ++i) { 2308 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 2309 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 2310 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 2311 } 2312 2313 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 2314 PetscObject obj; 2315 PetscClassId id; 2316 PetscReal *B; 2317 PetscInt NcJ = 0, cpdim = 0, j, qNc; 2318 2319 ierr = PetscDSGetDiscretization(prob, fieldJ, &obj);CHKERRQ(ierr); 2320 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2321 if (id == PETSCFE_CLASSID) { 2322 PetscFE fe = (PetscFE) obj; 2323 2324 /* Evaluate basis at points */ 2325 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 2326 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 2327 /* For now, fields only interpolate themselves */ 2328 if (fieldI == fieldJ) { 2329 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); 2330 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 2331 for (i = 0, k = 0; i < fpdim; ++i) { 2332 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 2333 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);CHKERRQ(ierr); 2334 if (qNc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in quadrature %D does not match coarse field %D", qNc, NcJ); 2335 for (p = 0; p < Np; ++p, ++k) { 2336 for (j = 0; j < cpdim; ++j) { 2337 /* 2338 cTotDim: Total columns in element interpolation matrix, sum of number of dual basis functionals in each field 2339 offsetI, offsetJ: Offsets into the larger element interpolation matrix for different fields 2340 fpdim, i, cpdim, j: Dofs for fine and coarse grids, correspond to dual space basis functionals 2341 qNC, Nc, Ncj, c: Number of components in this field 2342 Np, p: Number of quad points in the fine grid functional i 2343 k: i*Np + p, overall point number for the interpolation 2344 */ 2345 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p*qNc+c]; 2346 } 2347 } 2348 } 2349 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 2350 } 2351 } else if (id == PETSCFV_CLASSID) { 2352 PetscFV fv = (PetscFV) obj; 2353 2354 /* Evaluate constant function at points */ 2355 ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr); 2356 cpdim = 1; 2357 /* For now, fields only interpolate themselves */ 2358 if (fieldI == fieldJ) { 2359 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); 2360 for (i = 0, k = 0; i < fpdim; ++i) { 2361 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 2362 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);CHKERRQ(ierr); 2363 if (qNc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in quadrature %D does not match coarse field %D", qNc, NcJ); 2364 for (p = 0; p < Np; ++p, ++k) { 2365 for (j = 0; j < cpdim; ++j) { 2366 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += 1.0*qweights[p*qNc+c]; 2367 } 2368 } 2369 } 2370 } 2371 } 2372 offsetJ += cpdim; 2373 } 2374 offsetI += fpdim; 2375 ierr = PetscFree(points);CHKERRQ(ierr); 2376 } 2377 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 2378 /* Preallocate matrix */ 2379 { 2380 Mat preallocator; 2381 PetscScalar *vals; 2382 PetscInt *cellCIndices, *cellFIndices; 2383 PetscInt locRows, locCols, cell; 2384 2385 ierr = MatGetLocalSize(In, &locRows, &locCols);CHKERRQ(ierr); 2386 ierr = MatCreate(PetscObjectComm((PetscObject) In), &preallocator);CHKERRQ(ierr); 2387 ierr = MatSetType(preallocator, MATPREALLOCATOR);CHKERRQ(ierr); 2388 ierr = MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 2389 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 2390 ierr = PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 2391 for (cell = cStart; cell < cEnd; ++cell) { 2392 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 2393 ierr = MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);CHKERRQ(ierr); 2394 } 2395 ierr = PetscFree3(vals,cellCIndices,cellFIndices);CHKERRQ(ierr); 2396 ierr = MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2397 ierr = MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2398 ierr = MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);CHKERRQ(ierr); 2399 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 2400 } 2401 /* Fill matrix */ 2402 ierr = MatZeroEntries(In);CHKERRQ(ierr); 2403 for (c = cStart; c < cEnd; ++c) { 2404 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 2405 } 2406 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 2407 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 2408 ierr = PetscFree(elemMat);CHKERRQ(ierr); 2409 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2410 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2411 if (mesh->printFEM) { 2412 ierr = PetscPrintf(PetscObjectComm((PetscObject)In), "%s:\n", name);CHKERRQ(ierr); 2413 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 2414 ierr = MatView(In, NULL);CHKERRQ(ierr); 2415 } 2416 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2417 PetscFunctionReturn(0); 2418 } 2419 2420 PetscErrorCode DMPlexComputeMassMatrixNested(DM dmc, DM dmf, Mat mass, void *user) 2421 { 2422 SETERRQ(PetscObjectComm((PetscObject) dmc), PETSC_ERR_SUP, "Laziness"); 2423 } 2424 2425 /*@ 2426 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 2427 2428 Input Parameters: 2429 + dmf - The fine mesh 2430 . dmc - The coarse mesh 2431 - user - The user context 2432 2433 Output Parameter: 2434 . In - The interpolation matrix 2435 2436 Level: developer 2437 2438 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 2439 @*/ 2440 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 2441 { 2442 DM_Plex *mesh = (DM_Plex *) dmf->data; 2443 const char *name = "Interpolator"; 2444 PetscDS prob; 2445 PetscSection fsection, csection, globalFSection, globalCSection; 2446 PetscHSetIJ ht; 2447 PetscLayout rLayout; 2448 PetscInt *dnz, *onz; 2449 PetscInt locRows, rStart, rEnd; 2450 PetscReal *x, *v0, *J, *invJ, detJ; 2451 PetscReal *v0c, *Jc, *invJc, detJc; 2452 PetscScalar *elemMat; 2453 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 2454 PetscErrorCode ierr; 2455 2456 PetscFunctionBegin; 2457 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2458 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 2459 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 2460 ierr = PetscDSGetWorkspace(prob, &x, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 2461 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 2462 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 2463 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 2464 ierr = DMGetSection(dmf, &fsection);CHKERRQ(ierr); 2465 ierr = DMGetGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 2466 ierr = DMGetSection(dmc, &csection);CHKERRQ(ierr); 2467 ierr = DMGetGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 2468 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 2469 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 2470 ierr = PetscMalloc1(totDim, &elemMat);CHKERRQ(ierr); 2471 2472 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 2473 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 2474 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 2475 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 2476 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 2477 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 2478 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 2479 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 2480 ierr = PetscHSetIJCreate(&ht);CHKERRQ(ierr); 2481 for (field = 0; field < Nf; ++field) { 2482 PetscObject obj; 2483 PetscClassId id; 2484 PetscDualSpace Q = NULL; 2485 PetscQuadrature f; 2486 const PetscReal *qpoints; 2487 PetscInt Nc, Np, fpdim, i, d; 2488 2489 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 2490 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2491 if (id == PETSCFE_CLASSID) { 2492 PetscFE fe = (PetscFE) obj; 2493 2494 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 2495 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 2496 } else if (id == PETSCFV_CLASSID) { 2497 PetscFV fv = (PetscFV) obj; 2498 2499 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 2500 Nc = 1; 2501 } 2502 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 2503 /* For each fine grid cell */ 2504 for (cell = cStart; cell < cEnd; ++cell) { 2505 PetscInt *findices, *cindices; 2506 PetscInt numFIndices, numCIndices; 2507 2508 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2509 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 2510 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim); 2511 for (i = 0; i < fpdim; ++i) { 2512 Vec pointVec; 2513 PetscScalar *pV; 2514 PetscSF coarseCellSF = NULL; 2515 const PetscSFNode *coarseCells; 2516 PetscInt numCoarseCells, q, c; 2517 2518 /* Get points from the dual basis functional quadrature */ 2519 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 2520 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 2521 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 2522 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 2523 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2524 for (q = 0; q < Np; ++q) { 2525 const PetscReal xi0[3] = {-1., -1., -1.}; 2526 2527 /* Transform point to real space */ 2528 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x); 2529 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 2530 } 2531 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2532 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 2533 /* OPT: Pack all quad points from fine cell */ 2534 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 2535 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 2536 /* Update preallocation info */ 2537 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 2538 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 2539 { 2540 PetscHashIJKey key; 2541 PetscBool missing; 2542 2543 key.i = findices[i]; 2544 if (key.i >= 0) { 2545 /* Get indices for coarse elements */ 2546 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 2547 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2548 for (c = 0; c < numCIndices; ++c) { 2549 key.j = cindices[c]; 2550 if (key.j < 0) continue; 2551 ierr = PetscHSetIJQueryAdd(ht, key, &missing);CHKERRQ(ierr); 2552 if (missing) { 2553 if ((key.j >= rStart) && (key.j < rEnd)) ++dnz[key.i-rStart]; 2554 else ++onz[key.i-rStart]; 2555 } 2556 } 2557 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2558 } 2559 } 2560 } 2561 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 2562 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 2563 } 2564 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2565 } 2566 } 2567 ierr = PetscHSetIJDestroy(&ht);CHKERRQ(ierr); 2568 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 2569 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2570 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 2571 for (field = 0; field < Nf; ++field) { 2572 PetscObject obj; 2573 PetscClassId id; 2574 PetscDualSpace Q = NULL; 2575 PetscQuadrature f; 2576 const PetscReal *qpoints, *qweights; 2577 PetscInt Nc, qNc, Np, fpdim, i, d; 2578 2579 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 2580 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2581 if (id == PETSCFE_CLASSID) { 2582 PetscFE fe = (PetscFE) obj; 2583 2584 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 2585 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 2586 } else if (id == PETSCFV_CLASSID) { 2587 PetscFV fv = (PetscFV) obj; 2588 2589 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 2590 Nc = 1; 2591 } else SETERRQ1(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONG,"Unknown discretization type for field %d",field); 2592 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 2593 /* For each fine grid cell */ 2594 for (cell = cStart; cell < cEnd; ++cell) { 2595 PetscInt *findices, *cindices; 2596 PetscInt numFIndices, numCIndices; 2597 2598 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2599 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 2600 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim); 2601 for (i = 0; i < fpdim; ++i) { 2602 Vec pointVec; 2603 PetscScalar *pV; 2604 PetscSF coarseCellSF = NULL; 2605 const PetscSFNode *coarseCells; 2606 PetscInt numCoarseCells, cpdim, q, c, j; 2607 2608 /* Get points from the dual basis functional quadrature */ 2609 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 2610 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);CHKERRQ(ierr); 2611 if (qNc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in quadrature %D does not match coarse field %D", qNc, Nc); 2612 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 2613 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 2614 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2615 for (q = 0; q < Np; ++q) { 2616 const PetscReal xi0[3] = {-1., -1., -1.}; 2617 2618 /* Transform point to real space */ 2619 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x); 2620 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 2621 } 2622 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2623 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 2624 /* OPT: Read this out from preallocation information */ 2625 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 2626 /* Update preallocation info */ 2627 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 2628 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 2629 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2630 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 2631 PetscReal pVReal[3]; 2632 const PetscReal xi0[3] = {-1., -1., -1.}; 2633 2634 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2635 /* Transform points from real space to coarse reference space */ 2636 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 2637 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 2638 CoordinatesRealToRef(dim, dim, xi0, v0c, invJc, pVReal, x); 2639 2640 if (id == PETSCFE_CLASSID) { 2641 PetscFE fe = (PetscFE) obj; 2642 PetscReal *B; 2643 2644 /* Evaluate coarse basis on contained point */ 2645 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 2646 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 2647 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 2648 /* Get elemMat entries by multiplying by weight */ 2649 for (j = 0; j < cpdim; ++j) { 2650 for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*qweights[ccell*qNc + c]; 2651 } 2652 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 2653 } else { 2654 cpdim = 1; 2655 for (j = 0; j < cpdim; ++j) { 2656 for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c]; 2657 } 2658 } 2659 /* Update interpolator */ 2660 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 2661 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 2662 ierr = MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 2663 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2664 } 2665 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2666 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 2667 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 2668 } 2669 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2670 } 2671 } 2672 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 2673 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 2674 ierr = PetscFree(elemMat);CHKERRQ(ierr); 2675 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2676 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2677 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2678 PetscFunctionReturn(0); 2679 } 2680 2681 /*@ 2682 DMPlexComputeMassMatrixGeneral - Form the local portion of the mass matrix M from the coarse DM to a non-nested fine DM. 2683 2684 Input Parameters: 2685 + dmf - The fine mesh 2686 . dmc - The coarse mesh 2687 - user - The user context 2688 2689 Output Parameter: 2690 . mass - The mass matrix 2691 2692 Level: developer 2693 2694 .seealso: DMPlexComputeMassMatrixNested(), DMPlexComputeInterpolatorNested(), DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 2695 @*/ 2696 PetscErrorCode DMPlexComputeMassMatrixGeneral(DM dmc, DM dmf, Mat mass, void *user) 2697 { 2698 DM_Plex *mesh = (DM_Plex *) dmf->data; 2699 const char *name = "Mass Matrix"; 2700 PetscDS prob; 2701 PetscSection fsection, csection, globalFSection, globalCSection; 2702 PetscHSetIJ ht; 2703 PetscLayout rLayout; 2704 PetscInt *dnz, *onz; 2705 PetscInt locRows, rStart, rEnd; 2706 PetscReal *x, *v0, *J, *invJ, detJ; 2707 PetscReal *v0c, *Jc, *invJc, detJc; 2708 PetscScalar *elemMat; 2709 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 2710 PetscErrorCode ierr; 2711 2712 PetscFunctionBegin; 2713 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 2714 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 2715 ierr = PetscDSGetWorkspace(prob, &x, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 2716 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 2717 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 2718 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 2719 ierr = DMGetSection(dmf, &fsection);CHKERRQ(ierr); 2720 ierr = DMGetGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 2721 ierr = DMGetSection(dmc, &csection);CHKERRQ(ierr); 2722 ierr = DMGetGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 2723 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 2724 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 2725 ierr = PetscMalloc1(totDim, &elemMat);CHKERRQ(ierr); 2726 2727 ierr = MatGetLocalSize(mass, &locRows, NULL);CHKERRQ(ierr); 2728 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) mass), &rLayout);CHKERRQ(ierr); 2729 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 2730 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 2731 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 2732 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 2733 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 2734 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 2735 ierr = PetscHSetIJCreate(&ht);CHKERRQ(ierr); 2736 for (field = 0; field < Nf; ++field) { 2737 PetscObject obj; 2738 PetscClassId id; 2739 PetscQuadrature quad; 2740 const PetscReal *qpoints; 2741 PetscInt Nq, Nc, i, d; 2742 2743 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 2744 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2745 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetQuadrature((PetscFE) obj, &quad);CHKERRQ(ierr);} 2746 else {ierr = PetscFVGetQuadrature((PetscFV) obj, &quad);CHKERRQ(ierr);} 2747 ierr = PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, NULL);CHKERRQ(ierr); 2748 /* For each fine grid cell */ 2749 for (cell = cStart; cell < cEnd; ++cell) { 2750 Vec pointVec; 2751 PetscScalar *pV; 2752 PetscSF coarseCellSF = NULL; 2753 const PetscSFNode *coarseCells; 2754 PetscInt numCoarseCells, q, c; 2755 PetscInt *findices, *cindices; 2756 PetscInt numFIndices, numCIndices; 2757 2758 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2759 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 2760 /* Get points from the quadrature */ 2761 ierr = VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);CHKERRQ(ierr); 2762 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 2763 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2764 for (q = 0; q < Nq; ++q) { 2765 const PetscReal xi0[3] = {-1., -1., -1.}; 2766 2767 /* Transform point to real space */ 2768 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x); 2769 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 2770 } 2771 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2772 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 2773 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 2774 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 2775 /* Update preallocation info */ 2776 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 2777 if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 2778 { 2779 PetscHashIJKey key; 2780 PetscBool missing; 2781 2782 for (i = 0; i < numFIndices; ++i) { 2783 key.i = findices[i]; 2784 if (key.i >= 0) { 2785 /* Get indices for coarse elements */ 2786 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 2787 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2788 for (c = 0; c < numCIndices; ++c) { 2789 key.j = cindices[c]; 2790 if (key.j < 0) continue; 2791 ierr = PetscHSetIJQueryAdd(ht, key, &missing);CHKERRQ(ierr); 2792 if (missing) { 2793 if ((key.j >= rStart) && (key.j < rEnd)) ++dnz[key.i-rStart]; 2794 else ++onz[key.i-rStart]; 2795 } 2796 } 2797 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2798 } 2799 } 2800 } 2801 } 2802 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 2803 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 2804 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2805 } 2806 } 2807 ierr = PetscHSetIJDestroy(&ht);CHKERRQ(ierr); 2808 ierr = MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 2809 ierr = MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2810 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 2811 for (field = 0; field < Nf; ++field) { 2812 PetscObject obj; 2813 PetscClassId id; 2814 PetscQuadrature quad; 2815 PetscReal *Bfine; 2816 const PetscReal *qpoints, *qweights; 2817 PetscInt Nq, Nc, i, d; 2818 2819 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 2820 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2821 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetQuadrature((PetscFE) obj, &quad);CHKERRQ(ierr);ierr = PetscFEGetDefaultTabulation((PetscFE) obj, &Bfine, NULL, NULL);CHKERRQ(ierr);} 2822 else {ierr = PetscFVGetQuadrature((PetscFV) obj, &quad);CHKERRQ(ierr);} 2823 ierr = PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, &qweights);CHKERRQ(ierr); 2824 /* For each fine grid cell */ 2825 for (cell = cStart; cell < cEnd; ++cell) { 2826 Vec pointVec; 2827 PetscScalar *pV; 2828 PetscSF coarseCellSF = NULL; 2829 const PetscSFNode *coarseCells; 2830 PetscInt numCoarseCells, cpdim, q, c, j; 2831 PetscInt *findices, *cindices; 2832 PetscInt numFIndices, numCIndices; 2833 2834 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2835 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 2836 /* Get points from the quadrature */ 2837 ierr = VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);CHKERRQ(ierr); 2838 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 2839 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2840 for (q = 0; q < Nq; ++q) { 2841 const PetscReal xi0[3] = {-1., -1., -1.}; 2842 2843 /* Transform point to real space */ 2844 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x); 2845 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 2846 } 2847 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2848 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 2849 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 2850 /* Update matrix */ 2851 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 2852 if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 2853 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2854 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 2855 PetscReal pVReal[3]; 2856 const PetscReal xi0[3] = {-1., -1., -1.}; 2857 2858 2859 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2860 /* Transform points from real space to coarse reference space */ 2861 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 2862 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 2863 CoordinatesRealToRef(dim, dim, xi0, v0c, invJc, pVReal, x); 2864 2865 if (id == PETSCFE_CLASSID) { 2866 PetscFE fe = (PetscFE) obj; 2867 PetscReal *B; 2868 2869 /* Evaluate coarse basis on contained point */ 2870 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 2871 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 2872 /* Get elemMat entries by multiplying by weight */ 2873 for (i = 0; i < numFIndices; ++i) { 2874 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 2875 for (j = 0; j < cpdim; ++j) { 2876 for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*Bfine[(ccell*numFIndices + i)*Nc + c]*qweights[ccell*Nc + c]*detJ; 2877 } 2878 /* Update interpolator */ 2879 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 2880 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 2881 ierr = MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);CHKERRQ(ierr); 2882 } 2883 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 2884 } else { 2885 cpdim = 1; 2886 for (i = 0; i < numFIndices; ++i) { 2887 ierr = PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));CHKERRQ(ierr); 2888 for (j = 0; j < cpdim; ++j) { 2889 for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*1.0*qweights[ccell*Nc + c]*detJ; 2890 } 2891 /* Update interpolator */ 2892 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 2893 ierr = PetscPrintf(PETSC_COMM_SELF, "Nq: %d %d Nf: %d %d Nc: %d %d\n", ccell, Nq, i, numFIndices, j, numCIndices);CHKERRQ(ierr); 2894 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 2895 ierr = MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);CHKERRQ(ierr); 2896 } 2897 } 2898 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 2899 } 2900 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2901 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 2902 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 2903 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 2904 } 2905 } 2906 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 2907 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 2908 ierr = PetscFree(elemMat);CHKERRQ(ierr); 2909 ierr = MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2910 ierr = MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2911 PetscFunctionReturn(0); 2912 } 2913 2914 /*@ 2915 DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns 2916 2917 Input Parameters: 2918 + dmc - The coarse mesh 2919 - dmf - The fine mesh 2920 - user - The user context 2921 2922 Output Parameter: 2923 . sc - The mapping 2924 2925 Level: developer 2926 2927 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 2928 @*/ 2929 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 2930 { 2931 PetscDS prob; 2932 PetscFE *feRef; 2933 PetscFV *fvRef; 2934 Vec fv, cv; 2935 IS fis, cis; 2936 PetscSection fsection, fglobalSection, csection, cglobalSection; 2937 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 2938 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m; 2939 PetscBool *needAvg; 2940 PetscErrorCode ierr; 2941 2942 PetscFunctionBegin; 2943 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2944 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 2945 ierr = DMGetSection(dmf, &fsection);CHKERRQ(ierr); 2946 ierr = DMGetGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 2947 ierr = DMGetSection(dmc, &csection);CHKERRQ(ierr); 2948 ierr = DMGetGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 2949 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 2950 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 2951 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 2952 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 2953 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 2954 ierr = PetscCalloc3(Nf,&feRef,Nf,&fvRef,Nf,&needAvg);CHKERRQ(ierr); 2955 for (f = 0; f < Nf; ++f) { 2956 PetscObject obj; 2957 PetscClassId id; 2958 PetscInt fNb = 0, Nc = 0; 2959 2960 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 2961 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2962 if (id == PETSCFE_CLASSID) { 2963 PetscFE fe = (PetscFE) obj; 2964 PetscSpace sp; 2965 PetscInt maxDegree; 2966 2967 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 2968 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 2969 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 2970 ierr = PetscFEGetBasisSpace(fe, &sp);CHKERRQ(ierr); 2971 ierr = PetscSpaceGetDegree(sp, NULL, &maxDegree);CHKERRQ(ierr); 2972 if (!maxDegree) needAvg[f] = PETSC_TRUE; 2973 } else if (id == PETSCFV_CLASSID) { 2974 PetscFV fv = (PetscFV) obj; 2975 PetscDualSpace Q; 2976 2977 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 2978 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 2979 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 2980 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 2981 needAvg[f] = PETSC_TRUE; 2982 } 2983 fTotDim += fNb; 2984 } 2985 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 2986 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 2987 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 2988 PetscFE feC; 2989 PetscFV fvC; 2990 PetscDualSpace QF, QC; 2991 PetscInt order = -1, NcF, NcC, fpdim, cpdim; 2992 2993 if (feRef[field]) { 2994 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 2995 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 2996 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 2997 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 2998 ierr = PetscDualSpaceGetOrder(QF, &order);CHKERRQ(ierr); 2999 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 3000 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 3001 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 3002 } else { 3003 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 3004 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 3005 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 3006 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 3007 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 3008 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 3009 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 3010 } 3011 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); 3012 for (c = 0; c < cpdim; ++c) { 3013 PetscQuadrature cfunc; 3014 const PetscReal *cqpoints, *cqweights; 3015 PetscInt NqcC, NpC; 3016 PetscBool found = PETSC_FALSE; 3017 3018 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 3019 ierr = PetscQuadratureGetData(cfunc, NULL, &NqcC, &NpC, &cqpoints, &cqweights);CHKERRQ(ierr); 3020 if (NqcC != NcC) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of quadrature components %D must match number of field components", NqcC, NcC); 3021 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 3022 for (f = 0; f < fpdim; ++f) { 3023 PetscQuadrature ffunc; 3024 const PetscReal *fqpoints, *fqweights; 3025 PetscReal sum = 0.0; 3026 PetscInt NqcF, NpF; 3027 3028 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 3029 ierr = PetscQuadratureGetData(ffunc, NULL, &NqcF, &NpF, &fqpoints, &fqweights);CHKERRQ(ierr); 3030 if (NqcF != NcF) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of quadrature components %D must match number of field components", NqcF, NcF); 3031 if (NpC != NpF) continue; 3032 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 3033 if (sum > 1.0e-9) continue; 3034 for (d = 0; d < NcC; ++d) sum += PetscAbsReal(cqweights[d]*fqweights[d]); 3035 if (sum < 1.0e-9) continue; 3036 cmap[offsetC+c] = offsetF+f; 3037 found = PETSC_TRUE; 3038 break; 3039 } 3040 if (!found) { 3041 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 3042 if (fvRef[field] || (feRef[field] && order == 0)) { 3043 cmap[offsetC+c] = offsetF+0; 3044 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 3045 } 3046 } 3047 offsetC += cpdim; 3048 offsetF += fpdim; 3049 } 3050 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 3051 ierr = PetscFree3(feRef,fvRef,needAvg);CHKERRQ(ierr); 3052 3053 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 3054 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 3055 ierr = VecGetOwnershipRange(cv, &startC, &endC);CHKERRQ(ierr); 3056 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 3057 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 3058 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 3059 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 3060 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 3061 for (c = cStart; c < cEnd; ++c) { 3062 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 3063 for (d = 0; d < cTotDim; ++d) { 3064 if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue; 3065 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]]); 3066 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 3067 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 3068 } 3069 } 3070 ierr = PetscFree(cmap);CHKERRQ(ierr); 3071 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 3072 3073 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 3074 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 3075 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 3076 ierr = ISDestroy(&cis);CHKERRQ(ierr); 3077 ierr = ISDestroy(&fis);CHKERRQ(ierr); 3078 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 3079 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 3080 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 3081 PetscFunctionReturn(0); 3082 } 3083 3084 /*@C 3085 DMPlexGetCellFields - Retrieve the field values values for a chunk of cells 3086 3087 Input Parameters: 3088 + dm - The DM 3089 . cellIS - The cells to include 3090 . locX - A local vector with the solution fields 3091 . locX_t - A local vector with solution field time derivatives, or NULL 3092 - locA - A local vector with auxiliary fields, or NULL 3093 3094 Output Parameters: 3095 + u - The field coefficients 3096 . u_t - The fields derivative coefficients 3097 - a - The auxiliary field coefficients 3098 3099 Level: developer 3100 3101 .seealso: DMPlexGetFaceFields() 3102 @*/ 3103 PetscErrorCode DMPlexGetCellFields(DM dm, IS cellIS, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a) 3104 { 3105 DM plex, plexA = NULL; 3106 PetscSection section, sectionAux; 3107 PetscDS prob; 3108 const PetscInt *cells; 3109 PetscInt cStart, cEnd, numCells, totDim, totDimAux, c; 3110 PetscErrorCode ierr; 3111 3112 PetscFunctionBegin; 3113 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3114 PetscValidHeaderSpecific(locX, VEC_CLASSID, 4); 3115 if (locX_t) {PetscValidHeaderSpecific(locX_t, VEC_CLASSID, 5);} 3116 if (locA) {PetscValidHeaderSpecific(locA, VEC_CLASSID, 6);} 3117 PetscValidPointer(u, 7); 3118 PetscValidPointer(u_t, 8); 3119 PetscValidPointer(a, 9); 3120 ierr = DMPlexConvertPlex(dm, &plex, PETSC_FALSE);CHKERRQ(ierr); 3121 ierr = ISGetPointRange(cellIS, &cStart, &cEnd, &cells);CHKERRQ(ierr); 3122 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 3123 ierr = DMGetCellDS(dm, cStart, &prob);CHKERRQ(ierr); 3124 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 3125 if (locA) { 3126 DM dmAux; 3127 PetscDS probAux; 3128 3129 ierr = VecGetDM(locA, &dmAux);CHKERRQ(ierr); 3130 ierr = DMPlexConvertPlex(dmAux, &plexA, PETSC_FALSE);CHKERRQ(ierr); 3131 ierr = DMGetSection(dmAux, §ionAux);CHKERRQ(ierr); 3132 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 3133 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 3134 } 3135 numCells = cEnd - cStart; 3136 ierr = DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, u);CHKERRQ(ierr); 3137 if (locX_t) {ierr = DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, u_t);CHKERRQ(ierr);} else {*u_t = NULL;} 3138 if (locA) {ierr = DMGetWorkArray(dm, numCells*totDimAux, MPIU_SCALAR, a);CHKERRQ(ierr);} else {*a = NULL;} 3139 for (c = cStart; c < cEnd; ++c) { 3140 const PetscInt cell = cells ? cells[c] : c; 3141 const PetscInt cind = c - cStart; 3142 PetscScalar *x = NULL, *x_t = NULL, *ul = *u, *ul_t = *u_t, *al = *a; 3143 PetscInt i; 3144 3145 ierr = DMPlexVecGetClosure(plex, section, locX, cell, NULL, &x);CHKERRQ(ierr); 3146 for (i = 0; i < totDim; ++i) ul[cind*totDim+i] = x[i]; 3147 ierr = DMPlexVecRestoreClosure(plex, section, locX, cell, NULL, &x);CHKERRQ(ierr); 3148 if (locX_t) { 3149 ierr = DMPlexVecGetClosure(plex, section, locX_t, cell, NULL, &x_t);CHKERRQ(ierr); 3150 for (i = 0; i < totDim; ++i) ul_t[cind*totDim+i] = x_t[i]; 3151 ierr = DMPlexVecRestoreClosure(plex, section, locX_t, cell, NULL, &x_t);CHKERRQ(ierr); 3152 } 3153 if (locA) { 3154 PetscInt subcell; 3155 ierr = DMPlexGetAuxiliaryPoint(plex, plexA, cell, &subcell);CHKERRQ(ierr); 3156 ierr = DMPlexVecGetClosure(plexA, sectionAux, locA, subcell, NULL, &x);CHKERRQ(ierr); 3157 for (i = 0; i < totDimAux; ++i) al[cind*totDimAux+i] = x[i]; 3158 ierr = DMPlexVecRestoreClosure(plexA, sectionAux, locA, subcell, NULL, &x);CHKERRQ(ierr); 3159 } 3160 } 3161 ierr = DMDestroy(&plex);CHKERRQ(ierr); 3162 if (locA) {ierr = DMDestroy(&plexA);CHKERRQ(ierr);} 3163 ierr = ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);CHKERRQ(ierr); 3164 PetscFunctionReturn(0); 3165 } 3166 3167 /*@C 3168 DMPlexRestoreCellFields - Restore the field values values for a chunk of cells 3169 3170 Input Parameters: 3171 + dm - The DM 3172 . cellIS - The cells to include 3173 . locX - A local vector with the solution fields 3174 . locX_t - A local vector with solution field time derivatives, or NULL 3175 - locA - A local vector with auxiliary fields, or NULL 3176 3177 Output Parameters: 3178 + u - The field coefficients 3179 . u_t - The fields derivative coefficients 3180 - a - The auxiliary field coefficients 3181 3182 Level: developer 3183 3184 .seealso: DMPlexGetFaceFields() 3185 @*/ 3186 PetscErrorCode DMPlexRestoreCellFields(DM dm, IS cellIS, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a) 3187 { 3188 PetscErrorCode ierr; 3189 3190 PetscFunctionBegin; 3191 ierr = DMRestoreWorkArray(dm, 0, MPIU_SCALAR, u);CHKERRQ(ierr); 3192 if (locX_t) {ierr = DMRestoreWorkArray(dm, 0, MPIU_SCALAR, u_t);CHKERRQ(ierr);} 3193 if (locA) {ierr = DMRestoreWorkArray(dm, 0, MPIU_SCALAR, a);CHKERRQ(ierr);} 3194 PetscFunctionReturn(0); 3195 } 3196 3197 /*@C 3198 DMPlexGetFaceFields - Retrieve the field values values for a chunk of faces 3199 3200 Input Parameters: 3201 + dm - The DM 3202 . fStart - The first face to include 3203 . fEnd - The first face to exclude 3204 . locX - A local vector with the solution fields 3205 . locX_t - A local vector with solution field time derivatives, or NULL 3206 . faceGeometry - A local vector with face geometry 3207 . cellGeometry - A local vector with cell geometry 3208 - locaGrad - A local vector with field gradients, or NULL 3209 3210 Output Parameters: 3211 + Nface - The number of faces with field values 3212 . uL - The field values at the left side of the face 3213 - uR - The field values at the right side of the face 3214 3215 Level: developer 3216 3217 .seealso: DMPlexGetCellFields() 3218 @*/ 3219 PetscErrorCode DMPlexGetFaceFields(DM dm, PetscInt fStart, PetscInt fEnd, Vec locX, Vec locX_t, Vec faceGeometry, Vec cellGeometry, Vec locGrad, PetscInt *Nface, PetscScalar **uL, PetscScalar **uR) 3220 { 3221 DM dmFace, dmCell, dmGrad = NULL; 3222 PetscSection section; 3223 PetscDS prob; 3224 DMLabel ghostLabel; 3225 const PetscScalar *facegeom, *cellgeom, *x, *lgrad; 3226 PetscBool *isFE; 3227 PetscInt dim, Nf, f, Nc, numFaces = fEnd - fStart, iface, face; 3228 PetscErrorCode ierr; 3229 3230 PetscFunctionBegin; 3231 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3232 PetscValidHeaderSpecific(locX, VEC_CLASSID, 4); 3233 if (locX_t) {PetscValidHeaderSpecific(locX_t, VEC_CLASSID, 5);} 3234 PetscValidHeaderSpecific(faceGeometry, VEC_CLASSID, 6); 3235 PetscValidHeaderSpecific(cellGeometry, VEC_CLASSID, 7); 3236 if (locGrad) {PetscValidHeaderSpecific(locGrad, VEC_CLASSID, 8);} 3237 PetscValidPointer(uL, 9); 3238 PetscValidPointer(uR, 10); 3239 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3240 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 3241 ierr = DMGetSection(dm, §ion);CHKERRQ(ierr); 3242 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 3243 ierr = PetscDSGetTotalComponents(prob, &Nc);CHKERRQ(ierr); 3244 ierr = PetscMalloc1(Nf, &isFE);CHKERRQ(ierr); 3245 for (f = 0; f < Nf; ++f) { 3246 PetscObject obj; 3247 PetscClassId id; 3248 3249 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 3250 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 3251 if (id == PETSCFE_CLASSID) {isFE[f] = PETSC_TRUE;} 3252 else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;} 3253 else {isFE[f] = PETSC_FALSE;} 3254 } 3255 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 3256 ierr = VecGetArrayRead(locX, &x);CHKERRQ(ierr); 3257 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 3258 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 3259 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 3260 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 3261 if (locGrad) { 3262 ierr = VecGetDM(locGrad, &dmGrad);CHKERRQ(ierr); 3263 ierr = VecGetArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 3264 } 3265 ierr = DMGetWorkArray(dm, numFaces*Nc, MPIU_SCALAR, uL);CHKERRQ(ierr); 3266 ierr = DMGetWorkArray(dm, numFaces*Nc, MPIU_SCALAR, uR);CHKERRQ(ierr); 3267 /* Right now just eat the extra work for FE (could make a cell loop) */ 3268 for (face = fStart, iface = 0; face < fEnd; ++face) { 3269 const PetscInt *cells; 3270 PetscFVFaceGeom *fg; 3271 PetscFVCellGeom *cgL, *cgR; 3272 PetscScalar *xL, *xR, *gL, *gR; 3273 PetscScalar *uLl = *uL, *uRl = *uR; 3274 PetscInt ghost, nsupp, nchild; 3275 3276 ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr); 3277 ierr = DMPlexGetSupportSize(dm, face, &nsupp);CHKERRQ(ierr); 3278 ierr = DMPlexGetTreeChildren(dm, face, &nchild, NULL);CHKERRQ(ierr); 3279 if (ghost >= 0 || nsupp > 2 || nchild > 0) continue; 3280 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 3281 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 3282 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);CHKERRQ(ierr); 3283 ierr = DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);CHKERRQ(ierr); 3284 for (f = 0; f < Nf; ++f) { 3285 PetscInt off; 3286 3287 ierr = PetscDSGetComponentOffset(prob, f, &off);CHKERRQ(ierr); 3288 if (isFE[f]) { 3289 const PetscInt *cone; 3290 PetscInt comp, coneSizeL, coneSizeR, faceLocL, faceLocR, ldof, rdof, d; 3291 3292 xL = xR = NULL; 3293 ierr = PetscSectionGetFieldComponents(section, f, &comp);CHKERRQ(ierr); 3294 ierr = DMPlexVecGetClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);CHKERRQ(ierr); 3295 ierr = DMPlexVecGetClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);CHKERRQ(ierr); 3296 ierr = DMPlexGetCone(dm, cells[0], &cone);CHKERRQ(ierr); 3297 ierr = DMPlexGetConeSize(dm, cells[0], &coneSizeL);CHKERRQ(ierr); 3298 for (faceLocL = 0; faceLocL < coneSizeL; ++faceLocL) if (cone[faceLocL] == face) break; 3299 ierr = DMPlexGetCone(dm, cells[1], &cone);CHKERRQ(ierr); 3300 ierr = DMPlexGetConeSize(dm, cells[1], &coneSizeR);CHKERRQ(ierr); 3301 for (faceLocR = 0; faceLocR < coneSizeR; ++faceLocR) if (cone[faceLocR] == face) break; 3302 if (faceLocL == coneSizeL && faceLocR == coneSizeR) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %d in cone of cell %d or cell %d", face, cells[0], cells[1]); 3303 /* Check that FEM field has values in the right cell (sometimes its an FV ghost cell) */ 3304 /* TODO: this is a hack that might not be right for nonconforming */ 3305 if (faceLocL < coneSizeL) { 3306 ierr = PetscFEEvaluateFaceFields_Internal(prob, f, faceLocL, xL, &uLl[iface*Nc+off]);CHKERRQ(ierr); 3307 if (rdof == ldof && faceLocR < coneSizeR) {ierr = PetscFEEvaluateFaceFields_Internal(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);CHKERRQ(ierr);} 3308 else {for(d = 0; d < comp; ++d) uRl[iface*Nc+off+d] = uLl[iface*Nc+off+d];} 3309 } 3310 else { 3311 ierr = PetscFEEvaluateFaceFields_Internal(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);CHKERRQ(ierr); 3312 ierr = PetscSectionGetFieldComponents(section, f, &comp);CHKERRQ(ierr); 3313 for(d = 0; d < comp; ++d) uLl[iface*Nc+off+d] = uRl[iface*Nc+off+d]; 3314 } 3315 ierr = DMPlexVecRestoreClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);CHKERRQ(ierr); 3316 ierr = DMPlexVecRestoreClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);CHKERRQ(ierr); 3317 } else { 3318 PetscFV fv; 3319 PetscInt numComp, c; 3320 3321 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fv);CHKERRQ(ierr); 3322 ierr = PetscFVGetNumComponents(fv, &numComp);CHKERRQ(ierr); 3323 ierr = DMPlexPointLocalFieldRead(dm, cells[0], f, x, &xL);CHKERRQ(ierr); 3324 ierr = DMPlexPointLocalFieldRead(dm, cells[1], f, x, &xR);CHKERRQ(ierr); 3325 if (dmGrad) { 3326 PetscReal dxL[3], dxR[3]; 3327 3328 ierr = DMPlexPointLocalRead(dmGrad, cells[0], lgrad, &gL);CHKERRQ(ierr); 3329 ierr = DMPlexPointLocalRead(dmGrad, cells[1], lgrad, &gR);CHKERRQ(ierr); 3330 DMPlex_WaxpyD_Internal(dim, -1, cgL->centroid, fg->centroid, dxL); 3331 DMPlex_WaxpyD_Internal(dim, -1, cgR->centroid, fg->centroid, dxR); 3332 for (c = 0; c < numComp; ++c) { 3333 uLl[iface*Nc+off+c] = xL[c] + DMPlex_DotD_Internal(dim, &gL[c*dim], dxL); 3334 uRl[iface*Nc+off+c] = xR[c] + DMPlex_DotD_Internal(dim, &gR[c*dim], dxR); 3335 } 3336 } else { 3337 for (c = 0; c < numComp; ++c) { 3338 uLl[iface*Nc+off+c] = xL[c]; 3339 uRl[iface*Nc+off+c] = xR[c]; 3340 } 3341 } 3342 } 3343 } 3344 ++iface; 3345 } 3346 *Nface = iface; 3347 ierr = VecRestoreArrayRead(locX, &x);CHKERRQ(ierr); 3348 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 3349 ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 3350 if (locGrad) { 3351 ierr = VecRestoreArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 3352 } 3353 ierr = PetscFree(isFE);CHKERRQ(ierr); 3354 PetscFunctionReturn(0); 3355 } 3356 3357 /*@C 3358 DMPlexRestoreFaceFields - Restore the field values values for a chunk of faces 3359 3360 Input Parameters: 3361 + dm - The DM 3362 . fStart - The first face to include 3363 . fEnd - The first face to exclude 3364 . locX - A local vector with the solution fields 3365 . locX_t - A local vector with solution field time derivatives, or NULL 3366 . faceGeometry - A local vector with face geometry 3367 . cellGeometry - A local vector with cell geometry 3368 - locaGrad - A local vector with field gradients, or NULL 3369 3370 Output Parameters: 3371 + Nface - The number of faces with field values 3372 . uL - The field values at the left side of the face 3373 - uR - The field values at the right side of the face 3374 3375 Level: developer 3376 3377 .seealso: DMPlexGetFaceFields() 3378 @*/ 3379 PetscErrorCode DMPlexRestoreFaceFields(DM dm, PetscInt fStart, PetscInt fEnd, Vec locX, Vec locX_t, Vec faceGeometry, Vec cellGeometry, Vec locGrad, PetscInt *Nface, PetscScalar **uL, PetscScalar **uR) 3380 { 3381 PetscErrorCode ierr; 3382 3383 PetscFunctionBegin; 3384 ierr = DMRestoreWorkArray(dm, 0, MPIU_SCALAR, uL);CHKERRQ(ierr); 3385 ierr = DMRestoreWorkArray(dm, 0, MPIU_SCALAR, uR);CHKERRQ(ierr); 3386 PetscFunctionReturn(0); 3387 } 3388 3389 /*@C 3390 DMPlexGetFaceGeometry - Retrieve the geometric values for a chunk of faces 3391 3392 Input Parameters: 3393 + dm - The DM 3394 . fStart - The first face to include 3395 . fEnd - The first face to exclude 3396 . faceGeometry - A local vector with face geometry 3397 - cellGeometry - A local vector with cell geometry 3398 3399 Output Parameters: 3400 + Nface - The number of faces with field values 3401 . fgeom - The extract the face centroid and normal 3402 - vol - The cell volume 3403 3404 Level: developer 3405 3406 .seealso: DMPlexGetCellFields() 3407 @*/ 3408 PetscErrorCode DMPlexGetFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscInt *Nface, PetscFVFaceGeom **fgeom, PetscReal **vol) 3409 { 3410 DM dmFace, dmCell; 3411 DMLabel ghostLabel; 3412 const PetscScalar *facegeom, *cellgeom; 3413 PetscInt dim, numFaces = fEnd - fStart, iface, face; 3414 PetscErrorCode ierr; 3415 3416 PetscFunctionBegin; 3417 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3418 PetscValidHeaderSpecific(faceGeometry, VEC_CLASSID, 4); 3419 PetscValidHeaderSpecific(cellGeometry, VEC_CLASSID, 5); 3420 PetscValidPointer(fgeom, 6); 3421 PetscValidPointer(vol, 7); 3422 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3423 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 3424 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 3425 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 3426 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 3427 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 3428 ierr = PetscMalloc1(numFaces, fgeom);CHKERRQ(ierr); 3429 ierr = DMGetWorkArray(dm, numFaces*2, MPIU_SCALAR, vol);CHKERRQ(ierr); 3430 for (face = fStart, iface = 0; face < fEnd; ++face) { 3431 const PetscInt *cells; 3432 PetscFVFaceGeom *fg; 3433 PetscFVCellGeom *cgL, *cgR; 3434 PetscFVFaceGeom *fgeoml = *fgeom; 3435 PetscReal *voll = *vol; 3436 PetscInt ghost, d, nchild, nsupp; 3437 3438 ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr); 3439 ierr = DMPlexGetSupportSize(dm, face, &nsupp);CHKERRQ(ierr); 3440 ierr = DMPlexGetTreeChildren(dm, face, &nchild, NULL);CHKERRQ(ierr); 3441 if (ghost >= 0 || nsupp > 2 || nchild > 0) continue; 3442 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 3443 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 3444 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);CHKERRQ(ierr); 3445 ierr = DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);CHKERRQ(ierr); 3446 for (d = 0; d < dim; ++d) { 3447 fgeoml[iface].centroid[d] = fg->centroid[d]; 3448 fgeoml[iface].normal[d] = fg->normal[d]; 3449 } 3450 voll[iface*2+0] = cgL->volume; 3451 voll[iface*2+1] = cgR->volume; 3452 ++iface; 3453 } 3454 *Nface = iface; 3455 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 3456 ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 3457 PetscFunctionReturn(0); 3458 } 3459 3460 /*@C 3461 DMPlexRestoreFaceGeometry - Restore the field values values for a chunk of faces 3462 3463 Input Parameters: 3464 + dm - The DM 3465 . fStart - The first face to include 3466 . fEnd - The first face to exclude 3467 . faceGeometry - A local vector with face geometry 3468 - cellGeometry - A local vector with cell geometry 3469 3470 Output Parameters: 3471 + Nface - The number of faces with field values 3472 . fgeom - The extract the face centroid and normal 3473 - vol - The cell volume 3474 3475 Level: developer 3476 3477 .seealso: DMPlexGetFaceFields() 3478 @*/ 3479 PetscErrorCode DMPlexRestoreFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscInt *Nface, PetscFVFaceGeom **fgeom, PetscReal **vol) 3480 { 3481 PetscErrorCode ierr; 3482 3483 PetscFunctionBegin; 3484 ierr = PetscFree(*fgeom);CHKERRQ(ierr); 3485 ierr = DMRestoreWorkArray(dm, 0, MPIU_REAL, vol);CHKERRQ(ierr); 3486 PetscFunctionReturn(0); 3487 } 3488 3489 PetscErrorCode DMSNESGetFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom) 3490 { 3491 char composeStr[33] = {0}; 3492 PetscObjectId id; 3493 PetscContainer container; 3494 PetscErrorCode ierr; 3495 3496 PetscFunctionBegin; 3497 ierr = PetscObjectGetId((PetscObject)quad,&id);CHKERRQ(ierr); 3498 ierr = PetscSNPrintf(composeStr, 32, "DMSNESGetFEGeom_%x\n", id);CHKERRQ(ierr); 3499 ierr = PetscObjectQuery((PetscObject) pointIS, composeStr, (PetscObject *) &container);CHKERRQ(ierr); 3500 if (container) { 3501 ierr = PetscContainerGetPointer(container, (void **) geom);CHKERRQ(ierr); 3502 } else { 3503 ierr = DMFieldCreateFEGeom(coordField, pointIS, quad, faceData, geom);CHKERRQ(ierr); 3504 ierr = PetscContainerCreate(PETSC_COMM_SELF,&container);CHKERRQ(ierr); 3505 ierr = PetscContainerSetPointer(container, (void *) *geom);CHKERRQ(ierr); 3506 ierr = PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom);CHKERRQ(ierr); 3507 ierr = PetscObjectCompose((PetscObject) pointIS, composeStr, (PetscObject) container);CHKERRQ(ierr); 3508 ierr = PetscContainerDestroy(&container);CHKERRQ(ierr); 3509 } 3510 PetscFunctionReturn(0); 3511 } 3512 3513 PetscErrorCode DMSNESRestoreFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom) 3514 { 3515 PetscFunctionBegin; 3516 *geom = NULL; 3517 PetscFunctionReturn(0); 3518 } 3519 3520 PetscErrorCode DMPlexComputeResidual_Patch_Internal(DM dm, PetscSection section, IS cellIS, PetscReal t, Vec locX, Vec locX_t, Vec locF, void *user) 3521 { 3522 DM_Plex *mesh = (DM_Plex *) dm->data; 3523 const char *name = "Residual"; 3524 DM dmAux = NULL; 3525 DMLabel ghostLabel = NULL; 3526 PetscDS prob = NULL; 3527 PetscDS probAux = NULL; 3528 PetscBool useFEM = PETSC_FALSE; 3529 PetscBool isImplicit = (locX_t || t == PETSC_MIN_REAL) ? PETSC_TRUE : PETSC_FALSE; 3530 DMField coordField = NULL; 3531 Vec locA; 3532 PetscScalar *u = NULL, *u_t, *a, *uL = NULL, *uR = NULL; 3533 IS chunkIS; 3534 const PetscInt *cells; 3535 PetscInt cStart, cEnd, numCells; 3536 PetscInt Nf, f, totDim, totDimAux, numChunks, cellChunkSize, chunk, fStart, fEnd; 3537 PetscInt maxDegree = PETSC_MAX_INT; 3538 PetscQuadrature affineQuad = NULL, *quads = NULL; 3539 PetscFEGeom *affineGeom = NULL, **geoms = NULL; 3540 PetscErrorCode ierr; 3541 3542 PetscFunctionBegin; 3543 ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr); 3544 /* FEM+FVM */ 3545 /* 1: Get sizes from dm and dmAux */ 3546 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 3547 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 3548 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 3549 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 3550 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);CHKERRQ(ierr); 3551 if (locA) { 3552 ierr = VecGetDM(locA, &dmAux);CHKERRQ(ierr); 3553 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 3554 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 3555 } 3556 /* 2: Get geometric data */ 3557 for (f = 0; f < Nf; ++f) { 3558 PetscObject obj; 3559 PetscClassId id; 3560 PetscBool fimp; 3561 3562 ierr = PetscDSGetImplicit(prob, f, &fimp);CHKERRQ(ierr); 3563 if (isImplicit != fimp) continue; 3564 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 3565 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 3566 if (id == PETSCFE_CLASSID) {useFEM = PETSC_TRUE;} 3567 if (id == PETSCFV_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Use of FVM with PCPATCH not yet implemented"); 3568 } 3569 if (useFEM) { 3570 ierr = DMGetCoordinateField(dm, &coordField);CHKERRQ(ierr); 3571 ierr = DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);CHKERRQ(ierr); 3572 if (maxDegree <= 1) { 3573 ierr = DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);CHKERRQ(ierr); 3574 if (affineQuad) { 3575 ierr = DMSNESGetFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);CHKERRQ(ierr); 3576 } 3577 } else { 3578 ierr = PetscCalloc2(Nf,&quads,Nf,&geoms);CHKERRQ(ierr); 3579 for (f = 0; f < Nf; ++f) { 3580 PetscObject obj; 3581 PetscClassId id; 3582 PetscBool fimp; 3583 3584 ierr = PetscDSGetImplicit(prob, f, &fimp);CHKERRQ(ierr); 3585 if (isImplicit != fimp) continue; 3586 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 3587 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 3588 if (id == PETSCFE_CLASSID) { 3589 PetscFE fe = (PetscFE) obj; 3590 3591 ierr = PetscFEGetQuadrature(fe, &quads[f]);CHKERRQ(ierr); 3592 ierr = PetscObjectReference((PetscObject)quads[f]);CHKERRQ(ierr); 3593 ierr = DMSNESGetFEGeom(coordField,cellIS,quads[f],PETSC_FALSE,&geoms[f]);CHKERRQ(ierr); 3594 } 3595 } 3596 } 3597 } 3598 /* Loop over chunks */ 3599 ierr = ISGetPointRange(cellIS, &cStart, &cEnd, &cells);CHKERRQ(ierr); 3600 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 3601 if (useFEM) {ierr = ISCreate(PETSC_COMM_SELF, &chunkIS);CHKERRQ(ierr);} 3602 numCells = cEnd - cStart; 3603 numChunks = 1; 3604 cellChunkSize = numCells/numChunks; 3605 numChunks = PetscMin(1,numCells); 3606 for (chunk = 0; chunk < numChunks; ++chunk) { 3607 PetscScalar *elemVec, *fluxL = NULL, *fluxR = NULL; 3608 PetscReal *vol = NULL; 3609 PetscFVFaceGeom *fgeom = NULL; 3610 PetscInt cS = cStart+chunk*cellChunkSize, cE = PetscMin(cS+cellChunkSize, cEnd), numCells = cE - cS, c; 3611 PetscInt numFaces = 0; 3612 3613 /* Extract field coefficients */ 3614 if (useFEM) { 3615 ierr = ISGetPointSubrange(chunkIS, cS, cE, cells);CHKERRQ(ierr); 3616 ierr = DMPlexGetCellFields(dm, chunkIS, locX, locX_t, locA, &u, &u_t, &a);CHKERRQ(ierr); 3617 ierr = DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);CHKERRQ(ierr); 3618 ierr = PetscMemzero(elemVec, numCells*totDim * sizeof(PetscScalar));CHKERRQ(ierr); 3619 } 3620 /* TODO We will interlace both our field coefficients (u, u_t, uL, uR, etc.) and our output (elemVec, fL, fR). I think this works */ 3621 /* Loop over fields */ 3622 for (f = 0; f < Nf; ++f) { 3623 PetscObject obj; 3624 PetscClassId id; 3625 PetscBool fimp; 3626 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 3627 3628 ierr = PetscDSGetImplicit(prob, f, &fimp);CHKERRQ(ierr); 3629 if (isImplicit != fimp) continue; 3630 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 3631 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 3632 if (id == PETSCFE_CLASSID) { 3633 PetscFE fe = (PetscFE) obj; 3634 PetscFEGeom *geom = affineGeom ? affineGeom : geoms[f]; 3635 PetscFEGeom *chunkGeom = NULL; 3636 PetscQuadrature quad = affineQuad ? affineQuad : quads[f]; 3637 PetscInt Nq, Nb; 3638 3639 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 3640 ierr = PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 3641 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 3642 blockSize = Nb; 3643 batchSize = numBlocks * blockSize; 3644 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 3645 numChunks = numCells / (numBatches*batchSize); 3646 Ne = numChunks*numBatches*batchSize; 3647 Nr = numCells % (numBatches*batchSize); 3648 offset = numCells - Nr; 3649 /* Integrate FE residual to get elemVec (need fields at quadrature points) */ 3650 /* For FV, I think we use a P0 basis and the cell coefficients (for subdivided cells, we can tweak the basis tabulation to be the indicator function) */ 3651 ierr = PetscFEGeomGetChunk(geom,0,offset,&chunkGeom);CHKERRQ(ierr); 3652 ierr = PetscFEIntegrateResidual(prob, f, Ne, chunkGeom, u, u_t, probAux, a, t, elemVec);CHKERRQ(ierr); 3653 ierr = PetscFEGeomGetChunk(geom,offset,numCells,&chunkGeom);CHKERRQ(ierr); 3654 ierr = PetscFEIntegrateResidual(prob, f, Nr, chunkGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, &elemVec[offset*totDim]);CHKERRQ(ierr); 3655 ierr = PetscFEGeomRestoreChunk(geom,offset,numCells,&chunkGeom);CHKERRQ(ierr); 3656 } else if (id == PETSCFV_CLASSID) { 3657 PetscFV fv = (PetscFV) obj; 3658 3659 Ne = numFaces; 3660 /* Riemann solve over faces (need fields at face centroids) */ 3661 /* We need to evaluate FE fields at those coordinates */ 3662 ierr = PetscFVIntegrateRHSFunction(fv, prob, f, Ne, fgeom, vol, uL, uR, fluxL, fluxR);CHKERRQ(ierr); 3663 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 3664 } 3665 /* Loop over domain */ 3666 if (useFEM) { 3667 /* Add elemVec to locX */ 3668 for (c = cS; c < cE; ++c) { 3669 const PetscInt cell = cells ? cells[c] : c; 3670 const PetscInt cind = c - cStart; 3671 3672 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(cell, name, totDim, &elemVec[cind*totDim]);CHKERRQ(ierr);} 3673 if (ghostLabel) { 3674 PetscInt ghostVal; 3675 3676 ierr = DMLabelGetValue(ghostLabel,cell,&ghostVal);CHKERRQ(ierr); 3677 if (ghostVal > 0) continue; 3678 } 3679 ierr = DMPlexVecSetClosure(dm, section, locF, cell, &elemVec[cind*totDim], ADD_ALL_VALUES);CHKERRQ(ierr); 3680 } 3681 } 3682 /* Handle time derivative */ 3683 if (locX_t) { 3684 PetscScalar *x_t, *fa; 3685 3686 ierr = VecGetArray(locF, &fa);CHKERRQ(ierr); 3687 ierr = VecGetArray(locX_t, &x_t);CHKERRQ(ierr); 3688 for (f = 0; f < Nf; ++f) { 3689 PetscFV fv; 3690 PetscObject obj; 3691 PetscClassId id; 3692 PetscInt pdim, d; 3693 3694 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 3695 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 3696 if (id != PETSCFV_CLASSID) continue; 3697 fv = (PetscFV) obj; 3698 ierr = PetscFVGetNumComponents(fv, &pdim);CHKERRQ(ierr); 3699 for (c = cS; c < cE; ++c) { 3700 const PetscInt cell = cells ? cells[c] : c; 3701 PetscScalar *u_t, *r; 3702 3703 if (ghostLabel) { 3704 PetscInt ghostVal; 3705 3706 ierr = DMLabelGetValue(ghostLabel, cell, &ghostVal);CHKERRQ(ierr); 3707 if (ghostVal > 0) continue; 3708 } 3709 ierr = DMPlexPointLocalFieldRead(dm, cell, f, x_t, &u_t);CHKERRQ(ierr); 3710 ierr = DMPlexPointLocalFieldRef(dm, cell, f, fa, &r);CHKERRQ(ierr); 3711 for (d = 0; d < pdim; ++d) r[d] += u_t[d]; 3712 } 3713 } 3714 ierr = VecRestoreArray(locX_t, &x_t);CHKERRQ(ierr); 3715 ierr = VecRestoreArray(locF, &fa);CHKERRQ(ierr); 3716 } 3717 if (useFEM) { 3718 ierr = DMPlexRestoreCellFields(dm, chunkIS, locX, locX_t, locA, &u, &u_t, &a);CHKERRQ(ierr); 3719 ierr = DMRestoreWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);CHKERRQ(ierr); 3720 } 3721 } 3722 if (useFEM) {ierr = ISDestroy(&chunkIS);CHKERRQ(ierr);} 3723 ierr = ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);CHKERRQ(ierr); 3724 /* TODO Could include boundary residual here (see DMPlexComputeResidual_Internal) */ 3725 if (useFEM) { 3726 if (maxDegree <= 1) { 3727 ierr = DMSNESRestoreFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);CHKERRQ(ierr); 3728 ierr = PetscQuadratureDestroy(&affineQuad);CHKERRQ(ierr); 3729 } else { 3730 for (f = 0; f < Nf; ++f) { 3731 ierr = DMSNESRestoreFEGeom(coordField,cellIS,quads[f],PETSC_FALSE,&geoms[f]);CHKERRQ(ierr); 3732 ierr = PetscQuadratureDestroy(&quads[f]);CHKERRQ(ierr); 3733 } 3734 ierr = PetscFree2(quads,geoms);CHKERRQ(ierr); 3735 } 3736 } 3737 ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr); 3738 PetscFunctionReturn(0); 3739 } 3740 3741 /* 3742 We always assemble JacP, and if the matrix is different from Jac and two different sets of point functions are provided, we also assemble Jac 3743 3744 X - The local solution vector 3745 X_t - The local solution time derviative vector, or NULL 3746 */ 3747 PetscErrorCode DMPlexComputeJacobian_Patch_Internal(DM dm, PetscSection section, PetscSection globalSection, IS cellIS, 3748 PetscReal t, PetscReal X_tShift, Vec X, Vec X_t, Mat Jac, Mat JacP, void *ctx) 3749 { 3750 DM_Plex *mesh = (DM_Plex *) dm->data; 3751 const char *name = "Jacobian", *nameP = "JacobianPre"; 3752 DM dmAux = NULL; 3753 PetscDS prob, probAux = NULL; 3754 PetscSection sectionAux = NULL; 3755 Vec A; 3756 DMField coordField; 3757 PetscFEGeom *cgeomFEM; 3758 PetscQuadrature qGeom = NULL; 3759 Mat J = Jac, JP = JacP; 3760 PetscScalar *work, *u = NULL, *u_t = NULL, *a = NULL, *elemMat = NULL, *elemMatP = NULL, *elemMatD = NULL; 3761 PetscBool hasJac, hasPrec, hasDyn, assembleJac, isMatIS, isMatISP, *isFE, hasFV = PETSC_FALSE; 3762 const PetscInt *cells; 3763 PetscInt Nf, fieldI, fieldJ, maxDegree, numCells, cStart, cEnd, numChunks, chunkSize, chunk, totDim, totDimAux = 0, sz, wsz, off = 0, offCell = 0; 3764 PetscErrorCode ierr; 3765 3766 PetscFunctionBegin; 3767 CHKMEMQ; 3768 ierr = ISGetLocalSize(cellIS, &numCells);CHKERRQ(ierr); 3769 ierr = ISGetPointRange(cellIS, &cStart, &cEnd, &cells);CHKERRQ(ierr); 3770 ierr = PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr); 3771 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 3772 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 3773 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 3774 if (dmAux) { 3775 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 3776 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 3777 } 3778 /* Get flags */ 3779 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 3780 ierr = DMGetWorkArray(dm, Nf, MPIU_BOOL, &isFE);CHKERRQ(ierr); 3781 for (fieldI = 0; fieldI < Nf; ++fieldI) { 3782 PetscObject disc; 3783 PetscClassId id; 3784 ierr = PetscDSGetDiscretization(prob, fieldI, &disc);CHKERRQ(ierr); 3785 ierr = PetscObjectGetClassId(disc, &id);CHKERRQ(ierr); 3786 if (id == PETSCFE_CLASSID) {isFE[fieldI] = PETSC_TRUE;} 3787 else if (id == PETSCFV_CLASSID) {hasFV = PETSC_TRUE; isFE[fieldI] = PETSC_FALSE;} 3788 } 3789 ierr = PetscDSHasJacobian(prob, &hasJac);CHKERRQ(ierr); 3790 ierr = PetscDSHasJacobianPreconditioner(prob, &hasPrec);CHKERRQ(ierr); 3791 ierr = PetscDSHasDynamicJacobian(prob, &hasDyn);CHKERRQ(ierr); 3792 assembleJac = hasJac && hasPrec && (Jac != JacP) ? PETSC_TRUE : PETSC_FALSE; 3793 hasDyn = hasDyn && (X_tShift != 0.0) ? PETSC_TRUE : PETSC_FALSE; 3794 ierr = PetscObjectTypeCompare((PetscObject) Jac, MATIS, &isMatIS);CHKERRQ(ierr); 3795 ierr = PetscObjectTypeCompare((PetscObject) JacP, MATIS, &isMatISP);CHKERRQ(ierr); 3796 /* Setup input data and temp arrays (should be DMGetWorkArray) */ 3797 if (isMatISP || isMatISP) {ierr = DMPlexGetSubdomainSection(dm, &globalSection);CHKERRQ(ierr);} 3798 if (isMatIS) {ierr = MatISGetLocalMat(Jac, &J);CHKERRQ(ierr);} 3799 if (isMatISP) {ierr = MatISGetLocalMat(JacP, &JP);CHKERRQ(ierr);} 3800 if (hasFV) {ierr = MatSetOption(JP, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);CHKERRQ(ierr);} /* No allocated space for FV stuff, so ignore the zero entries */ 3801 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 3802 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 3803 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 3804 if (probAux) {ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr);} 3805 CHKMEMQ; 3806 /* Compute batch sizes */ 3807 if (isFE[0]) { 3808 PetscFE fe; 3809 PetscQuadrature q; 3810 PetscInt numQuadPoints, numBatches, batchSize, numBlocks, blockSize, Nb; 3811 3812 ierr = PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);CHKERRQ(ierr); 3813 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 3814 ierr = PetscQuadratureGetData(q, NULL, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 3815 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 3816 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 3817 blockSize = Nb*numQuadPoints; 3818 batchSize = numBlocks * blockSize; 3819 chunkSize = numBatches * batchSize; 3820 numChunks = numCells / chunkSize + numCells % chunkSize; 3821 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 3822 } else { 3823 chunkSize = numCells; 3824 numChunks = 1; 3825 } 3826 /* Get work space */ 3827 wsz = (((X?1:0) + (X_t?1:0) + (dmAux?1:0))*totDim + ((hasJac?1:0) + (hasPrec?1:0) + (hasDyn?1:0))*totDim*totDim)*chunkSize; 3828 ierr = DMGetWorkArray(dm, wsz, MPIU_SCALAR, &work);CHKERRQ(ierr); 3829 ierr = PetscMemzero(work, wsz * sizeof(PetscScalar));CHKERRQ(ierr); 3830 off = 0; 3831 u = X ? (sz = chunkSize*totDim, off += sz, work+off-sz) : NULL; 3832 u_t = X_t ? (sz = chunkSize*totDim, off += sz, work+off-sz) : NULL; 3833 a = dmAux ? (sz = chunkSize*totDimAux, off += sz, work+off-sz) : NULL; 3834 elemMat = hasJac ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL; 3835 elemMatP = hasPrec ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL; 3836 elemMatD = hasDyn ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL; 3837 if (off != wsz) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error is workspace size %D should be %D", off, wsz); 3838 /* Setup geometry */ 3839 ierr = DMGetCoordinateField(dm, &coordField);CHKERRQ(ierr); 3840 ierr = DMFieldGetDegree(coordField, cellIS, NULL, &maxDegree);CHKERRQ(ierr); 3841 if (maxDegree <= 1) {ierr = DMFieldCreateDefaultQuadrature(coordField, cellIS, &qGeom);CHKERRQ(ierr);} 3842 if (!qGeom) { 3843 PetscFE fe; 3844 3845 ierr = PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);CHKERRQ(ierr); 3846 ierr = PetscFEGetQuadrature(fe, &qGeom);CHKERRQ(ierr); 3847 ierr = PetscObjectReference((PetscObject) qGeom);CHKERRQ(ierr); 3848 } 3849 ierr = DMSNESGetFEGeom(coordField, cellIS, qGeom, PETSC_FALSE, &cgeomFEM);CHKERRQ(ierr); 3850 /* Compute volume integrals */ 3851 if (assembleJac) {ierr = MatZeroEntries(J);CHKERRQ(ierr);} 3852 ierr = MatZeroEntries(JP);CHKERRQ(ierr); 3853 for (chunk = 0; chunk < numChunks; ++chunk, offCell += chunkSize) { 3854 const PetscInt Ncell = PetscMin(chunkSize, numCells - offCell); 3855 PetscInt c; 3856 3857 /* Extract values */ 3858 for (c = 0; c < Ncell; ++c) { 3859 const PetscInt cell = cells ? cells[c+offCell] : c+offCell; 3860 PetscScalar *x = NULL, *x_t = NULL; 3861 PetscInt i; 3862 3863 if (X) { 3864 ierr = DMPlexVecGetClosure(dm, section, X, cell, NULL, &x);CHKERRQ(ierr); 3865 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 3866 ierr = DMPlexVecRestoreClosure(dm, section, X, cell, NULL, &x);CHKERRQ(ierr); 3867 } 3868 if (X_t) { 3869 ierr = DMPlexVecGetClosure(dm, section, X_t, cell, NULL, &x_t);CHKERRQ(ierr); 3870 for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i]; 3871 ierr = DMPlexVecRestoreClosure(dm, section, X_t, cell, NULL, &x_t);CHKERRQ(ierr); 3872 } 3873 if (dmAux) { 3874 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, cell, NULL, &x);CHKERRQ(ierr); 3875 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 3876 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, cell, NULL, &x);CHKERRQ(ierr); 3877 } 3878 } 3879 CHKMEMQ; 3880 for (fieldI = 0; fieldI < Nf; ++fieldI) { 3881 PetscFE fe; 3882 ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr); 3883 for (fieldJ = 0; fieldJ < Nf; ++fieldJ) { 3884 if (hasJac) {ierr = PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMat);CHKERRQ(ierr);} 3885 if (hasPrec) {ierr = PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_PRE, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMatP);CHKERRQ(ierr);} 3886 if (hasDyn) {ierr = PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMatD);CHKERRQ(ierr);} 3887 } 3888 /* For finite volume, add the identity */ 3889 if (!isFE[fieldI]) { 3890 PetscFV fv; 3891 PetscInt eOffset = 0, Nc, fc, foff; 3892 3893 ierr = PetscDSGetFieldOffset(prob, fieldI, &foff);CHKERRQ(ierr); 3894 ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fv);CHKERRQ(ierr); 3895 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 3896 for (c = 0; c < chunkSize; ++c, eOffset += totDim*totDim) { 3897 for (fc = 0; fc < Nc; ++fc) { 3898 const PetscInt i = foff + fc; 3899 if (hasJac) {elemMat [eOffset+i*totDim+i] = 1.0;} 3900 if (hasPrec) {elemMatP[eOffset+i*totDim+i] = 1.0;} 3901 } 3902 } 3903 } 3904 } 3905 CHKMEMQ; 3906 /* Add contribution from X_t */ 3907 if (hasDyn) {for (c = 0; c < chunkSize*totDim*totDim; ++c) elemMat[c] += X_tShift*elemMatD[c];} 3908 /* Insert values into matrix */ 3909 for (c = 0; c < Ncell; ++c) { 3910 const PetscInt cell = cells ? cells[c+offCell] : c+offCell; 3911 if (mesh->printFEM > 1) { 3912 if (hasJac) {ierr = DMPrintCellMatrix(cell, name, totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);CHKERRQ(ierr);} 3913 if (hasPrec) {ierr = DMPrintCellMatrix(cell, nameP, totDim, totDim, &elemMatP[(c-cStart)*totDim*totDim]);CHKERRQ(ierr);} 3914 } 3915 if (assembleJac) {ierr = DMPlexMatSetClosure(dm, section, globalSection, Jac, cell, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);CHKERRQ(ierr);} 3916 ierr = DMPlexMatSetClosure(dm, section, globalSection, JP, cell, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);CHKERRQ(ierr); 3917 } 3918 CHKMEMQ; 3919 } 3920 /* Cleanup */ 3921 ierr = DMSNESRestoreFEGeom(coordField, cellIS, qGeom, PETSC_FALSE, &cgeomFEM);CHKERRQ(ierr); 3922 ierr = PetscQuadratureDestroy(&qGeom);CHKERRQ(ierr); 3923 if (hasFV) {ierr = MatSetOption(JacP, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);CHKERRQ(ierr);} 3924 ierr = DMRestoreWorkArray(dm, Nf, MPIU_BOOL, &isFE);CHKERRQ(ierr); 3925 ierr = DMRestoreWorkArray(dm, ((1 + (X_t?1:0) + (dmAux?1:0))*totDim + ((hasJac?1:0) + (hasPrec?1:0) + (hasDyn?1:0))*totDim*totDim)*chunkSize, MPIU_SCALAR, &work);CHKERRQ(ierr); 3926 /* Compute boundary integrals */ 3927 /* ierr = DMPlexComputeBdJacobian_Internal(dm, X, X_t, t, X_tShift, Jac, JacP, ctx);CHKERRQ(ierr); */ 3928 /* Assemble matrix */ 3929 if (assembleJac) {ierr = MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);ierr = MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);} 3930 ierr = MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);ierr = MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3931 ierr = PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr); 3932 CHKMEMQ; 3933 PetscFunctionReturn(0); 3934 } 3935