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