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