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