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