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 PetscSection section; 1348 PetscQuadrature quad; 1349 Vec localX, tv; 1350 PetscFEGeom fegeom; 1351 PetscScalar *funcVal, *interpolant; 1352 PetscReal *coords, *gcoords; 1353 PetscReal *localDiff; 1354 const PetscReal *quadPoints, *quadWeights; 1355 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, c, field, fieldOffset; 1356 PetscBool transform; 1357 PetscErrorCode ierr; 1358 1359 PetscFunctionBegin; 1360 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1361 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 1362 ierr = DMGetLocalSection(dm, §ion);CHKERRQ(ierr); 1363 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1364 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1365 ierr = VecSet(localX, 0.0);CHKERRQ(ierr); 1366 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1367 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1368 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 1369 ierr = DMGetBasisTransformDM_Internal(dm, &tdm);CHKERRQ(ierr); 1370 ierr = DMGetBasisTransformVec_Internal(dm, &tv);CHKERRQ(ierr); 1371 ierr = DMHasBasisTransform(dm, &transform);CHKERRQ(ierr); 1372 for (field = 0; field < numFields; ++field) { 1373 PetscObject obj; 1374 PetscClassId id; 1375 PetscInt Nc; 1376 1377 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 1378 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1379 if (id == PETSCFE_CLASSID) { 1380 PetscFE fe = (PetscFE) obj; 1381 1382 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1383 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1384 } else if (id == PETSCFV_CLASSID) { 1385 PetscFV fv = (PetscFV) obj; 1386 1387 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1388 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1389 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field); 1390 numComponents += Nc; 1391 } 1392 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 1393 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 1394 ierr = PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,coordDim*(Nq+1),&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);CHKERRQ(ierr); 1395 ierr = DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1396 for (c = cStart; c < cEnd; ++c) { 1397 PetscScalar *x = NULL; 1398 PetscInt qc = 0; 1399 1400 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);CHKERRQ(ierr); 1401 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1402 1403 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1404 PetscObject obj; 1405 PetscClassId id; 1406 void * const ctx = ctxs ? ctxs[field] : NULL; 1407 PetscInt Nb, Nc, q, fc; 1408 1409 PetscReal elemDiff = 0.0; 1410 1411 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 1412 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1413 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1414 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1415 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field); 1416 if (debug) { 1417 char title[1024]; 1418 ierr = PetscSNPrintf(title, 1023, "Solution for Field %D", field);CHKERRQ(ierr); 1419 ierr = DMPrintCellVector(c, title, Nb, &x[fieldOffset]);CHKERRQ(ierr); 1420 } 1421 for (q = 0; q < Nq; ++q) { 1422 PetscFEGeom qgeom; 1423 1424 qgeom.dimEmbed = fegeom.dimEmbed; 1425 qgeom.J = &fegeom.J[q*coordDim*coordDim]; 1426 qgeom.invJ = &fegeom.invJ[q*coordDim*coordDim]; 1427 qgeom.detJ = &fegeom.detJ[q]; 1428 if (fegeom.detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature point %D", (double)fegeom.detJ[q], c, q); 1429 if (transform) { 1430 gcoords = &coords[coordDim*Nq]; 1431 ierr = DMPlexBasisTransformApplyReal_Internal(dm, &coords[coordDim*q], PETSC_TRUE, coordDim, &coords[coordDim*q], gcoords, dm->transformCtx);CHKERRQ(ierr); 1432 } else { 1433 gcoords = &coords[coordDim*q]; 1434 } 1435 ierr = (*funcs[field])(coordDim, time, gcoords, Nc, funcVal, ctx); 1436 if (ierr) { 1437 PetscErrorCode ierr2; 1438 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 1439 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 1440 ierr2 = PetscFree7(localDiff,funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2); 1441 CHKERRQ(ierr); 1442 } 1443 if (transform) {ierr = DMPlexBasisTransformApply_Internal(dm, &coords[coordDim*q], PETSC_FALSE, Nc, funcVal, funcVal, dm->transformCtx);CHKERRQ(ierr);} 1444 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], &qgeom, q, interpolant);CHKERRQ(ierr);} 1445 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1446 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field); 1447 for (fc = 0; fc < Nc; ++fc) { 1448 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 1449 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);} 1450 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q]; 1451 } 1452 } 1453 fieldOffset += Nb; 1454 qc += Nc; 1455 localDiff[field] += elemDiff; 1456 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %D field %D cum diff %g\n", c, field, (double)localDiff[field]);CHKERRQ(ierr);} 1457 } 1458 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1459 } 1460 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1461 ierr = MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1462 for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]); 1463 ierr = PetscFree7(localDiff,funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr); 1464 PetscFunctionReturn(0); 1465 } 1466 1467 /*@C 1468 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. 1469 1470 Collective on dm 1471 1472 Input Parameters: 1473 + dm - The DM 1474 . time - The time 1475 . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation 1476 . ctxs - Optional array of contexts to pass to each function, or NULL. 1477 - X - The coefficient vector u_h 1478 1479 Output Parameter: 1480 . D - A Vec which holds the difference ||u - u_h||_2 for each cell 1481 1482 Level: developer 1483 1484 .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff() 1485 @*/ 1486 PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D) 1487 { 1488 PetscSection section; 1489 PetscQuadrature quad; 1490 Vec localX; 1491 PetscFEGeom fegeom; 1492 PetscScalar *funcVal, *interpolant; 1493 PetscReal *coords; 1494 const PetscReal *quadPoints, *quadWeights; 1495 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, c, field, fieldOffset; 1496 PetscErrorCode ierr; 1497 1498 PetscFunctionBegin; 1499 ierr = VecSet(D, 0.0);CHKERRQ(ierr); 1500 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1501 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 1502 ierr = DMGetLocalSection(dm, §ion);CHKERRQ(ierr); 1503 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1504 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1505 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 1506 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1507 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1508 for (field = 0; field < numFields; ++field) { 1509 PetscObject obj; 1510 PetscClassId id; 1511 PetscInt Nc; 1512 1513 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 1514 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1515 if (id == PETSCFE_CLASSID) { 1516 PetscFE fe = (PetscFE) obj; 1517 1518 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1519 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1520 } else if (id == PETSCFV_CLASSID) { 1521 PetscFV fv = (PetscFV) obj; 1522 1523 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1524 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1525 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field); 1526 numComponents += Nc; 1527 } 1528 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 1529 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 1530 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);CHKERRQ(ierr); 1531 ierr = DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1532 for (c = cStart; c < cEnd; ++c) { 1533 PetscScalar *x = NULL; 1534 PetscScalar elemDiff = 0.0; 1535 PetscInt qc = 0; 1536 1537 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);CHKERRQ(ierr); 1538 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1539 1540 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1541 PetscObject obj; 1542 PetscClassId id; 1543 void * const ctx = ctxs ? ctxs[field] : NULL; 1544 PetscInt Nb, Nc, q, fc; 1545 1546 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 1547 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1548 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1549 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1550 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field); 1551 if (funcs[field]) { 1552 for (q = 0; q < Nq; ++q) { 1553 PetscFEGeom qgeom; 1554 1555 qgeom.dimEmbed = fegeom.dimEmbed; 1556 qgeom.J = &fegeom.J[q*coordDim*coordDim]; 1557 qgeom.invJ = &fegeom.invJ[q*coordDim*coordDim]; 1558 qgeom.detJ = &fegeom.detJ[q]; 1559 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); 1560 ierr = (*funcs[field])(coordDim, time, &coords[q*coordDim], Nc, funcVal, ctx); 1561 if (ierr) { 1562 PetscErrorCode ierr2; 1563 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 1564 ierr2 = PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2); 1565 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 1566 CHKERRQ(ierr); 1567 } 1568 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], &qgeom, q, interpolant);CHKERRQ(ierr);} 1569 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 1570 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field); 1571 for (fc = 0; fc < Nc; ++fc) { 1572 const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)]; 1573 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q]; 1574 } 1575 } 1576 } 1577 fieldOffset += Nb; 1578 qc += Nc; 1579 } 1580 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1581 ierr = VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);CHKERRQ(ierr); 1582 } 1583 ierr = PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr); 1584 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1585 ierr = VecSqrtAbs(D);CHKERRQ(ierr); 1586 PetscFunctionReturn(0); 1587 } 1588 1589 /*@C 1590 DMPlexComputeGradientClementInterpolant - This function computes the L2 projection of the cellwise gradient of a function u onto P1, and stores it in a Vec. 1591 1592 Collective on dm 1593 1594 Input Parameters: 1595 + dm - The DM 1596 - LocX - The coefficient vector u_h 1597 1598 Output Parameter: 1599 . locC - A Vec which holds the Clement interpolant of the gradient 1600 1601 Notes: 1602 Add citation to (Clement, 1975) and definition of the interpolant 1603 \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 1604 1605 Level: developer 1606 1607 .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff() 1608 @*/ 1609 PetscErrorCode DMPlexComputeGradientClementInterpolant(DM dm, Vec locX, Vec locC) 1610 { 1611 DM_Plex *mesh = (DM_Plex *) dm->data; 1612 PetscInt debug = mesh->printFEM; 1613 DM dmC; 1614 PetscSection section; 1615 PetscQuadrature quad; 1616 PetscScalar *interpolant, *gradsum; 1617 PetscFEGeom fegeom; 1618 PetscReal *coords; 1619 const PetscReal *quadPoints, *quadWeights; 1620 PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, vStart, vEnd, v, field, fieldOffset; 1621 PetscErrorCode ierr; 1622 1623 PetscFunctionBegin; 1624 ierr = VecGetDM(locC, &dmC);CHKERRQ(ierr); 1625 ierr = VecSet(locC, 0.0);CHKERRQ(ierr); 1626 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1627 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 1628 fegeom.dimEmbed = coordDim; 1629 ierr = DMGetLocalSection(dm, §ion);CHKERRQ(ierr); 1630 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1631 for (field = 0; field < numFields; ++field) { 1632 PetscObject obj; 1633 PetscClassId id; 1634 PetscInt Nc; 1635 1636 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 1637 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1638 if (id == PETSCFE_CLASSID) { 1639 PetscFE fe = (PetscFE) obj; 1640 1641 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1642 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1643 } else if (id == PETSCFV_CLASSID) { 1644 PetscFV fv = (PetscFV) obj; 1645 1646 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 1647 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1648 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field); 1649 numComponents += Nc; 1650 } 1651 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 1652 if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents); 1653 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); 1654 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1655 ierr = DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1656 for (v = vStart; v < vEnd; ++v) { 1657 PetscScalar volsum = 0.0; 1658 PetscInt *star = NULL; 1659 PetscInt starSize, st, d, fc; 1660 1661 ierr = PetscArrayzero(gradsum, coordDim*numComponents);CHKERRQ(ierr); 1662 ierr = DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1663 for (st = 0; st < starSize*2; st += 2) { 1664 const PetscInt cell = star[st]; 1665 PetscScalar *grad = &gradsum[coordDim*numComponents]; 1666 PetscScalar *x = NULL; 1667 PetscReal vol = 0.0; 1668 1669 if ((cell < cStart) || (cell >= cEnd)) continue; 1670 ierr = DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);CHKERRQ(ierr); 1671 ierr = DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr); 1672 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 1673 PetscObject obj; 1674 PetscClassId id; 1675 PetscInt Nb, Nc, q, qc = 0; 1676 1677 ierr = PetscArrayzero(grad, coordDim*numComponents);CHKERRQ(ierr); 1678 ierr = DMGetField(dm, field, NULL, &obj);CHKERRQ(ierr); 1679 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1680 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 1681 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 1682 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field); 1683 for (q = 0; q < Nq; ++q) { 1684 PetscFEGeom qgeom; 1685 1686 qgeom.dimEmbed = fegeom.dimEmbed; 1687 qgeom.J = &fegeom.J[q*coordDim*coordDim]; 1688 qgeom.invJ = &fegeom.invJ[q*coordDim*coordDim]; 1689 qgeom.detJ = &fegeom.detJ[q]; 1690 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); 1691 if (ierr) { 1692 PetscErrorCode ierr2; 1693 ierr2 = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr2); 1694 ierr2 = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr2); 1695 ierr2 = PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2); 1696 CHKERRQ(ierr); 1697 } 1698 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolateGradient_Static((PetscFE) obj, &x[fieldOffset], &qgeom, q, interpolant);CHKERRQ(ierr);} 1699 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field); 1700 for (fc = 0; fc < Nc; ++fc) { 1701 const PetscReal wt = quadWeights[q*qNc+qc+fc]; 1702 1703 for (d = 0; d < coordDim; ++d) grad[fc*coordDim+d] += interpolant[fc*dim+d]*wt*fegeom.detJ[q]; 1704 } 1705 vol += quadWeights[q*qNc]*fegeom.detJ[q]; 1706 } 1707 fieldOffset += Nb; 1708 qc += Nc; 1709 } 1710 ierr = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr); 1711 for (fc = 0; fc < numComponents; ++fc) { 1712 for (d = 0; d < coordDim; ++d) { 1713 gradsum[fc*coordDim+d] += grad[fc*coordDim+d]; 1714 } 1715 } 1716 volsum += vol; 1717 if (debug) { 1718 ierr = PetscPrintf(PETSC_COMM_SELF, "Cell %D gradient: [", cell);CHKERRQ(ierr); 1719 for (fc = 0; fc < numComponents; ++fc) { 1720 for (d = 0; d < coordDim; ++d) { 1721 if (fc || d > 0) {ierr = PetscPrintf(PETSC_COMM_SELF, ", ");CHKERRQ(ierr);} 1722 ierr = PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc*coordDim+d]));CHKERRQ(ierr); 1723 } 1724 } 1725 ierr = PetscPrintf(PETSC_COMM_SELF, "]\n");CHKERRQ(ierr); 1726 } 1727 } 1728 for (fc = 0; fc < numComponents; ++fc) { 1729 for (d = 0; d < coordDim; ++d) gradsum[fc*coordDim+d] /= volsum; 1730 } 1731 ierr = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 1732 ierr = DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES);CHKERRQ(ierr); 1733 } 1734 ierr = PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr); 1735 PetscFunctionReturn(0); 1736 } 1737 1738 static PetscErrorCode DMPlexComputeIntegral_Internal(DM dm, Vec X, PetscInt cStart, PetscInt cEnd, PetscScalar *cintegral, void *user) 1739 { 1740 DM dmAux = NULL; 1741 PetscDS prob, probAux = NULL; 1742 PetscSection section, sectionAux; 1743 Vec locX, locA; 1744 PetscInt dim, numCells = cEnd - cStart, c, f; 1745 PetscBool useFVM = PETSC_FALSE; 1746 /* DS */ 1747 PetscInt Nf, totDim, *uOff, *uOff_x, numConstants; 1748 PetscInt NfAux, totDimAux, *aOff; 1749 PetscScalar *u, *a; 1750 const PetscScalar *constants; 1751 /* Geometry */ 1752 PetscFEGeom *cgeomFEM; 1753 DM dmGrad; 1754 PetscQuadrature affineQuad = NULL; 1755 Vec cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL; 1756 PetscFVCellGeom *cgeomFVM; 1757 const PetscScalar *lgrad; 1758 PetscInt maxDegree; 1759 DMField coordField; 1760 IS cellIS; 1761 PetscErrorCode ierr; 1762 1763 PetscFunctionBegin; 1764 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1765 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1766 ierr = DMGetLocalSection(dm, §ion);CHKERRQ(ierr); 1767 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1768 /* Determine which discretizations we have */ 1769 for (f = 0; f < Nf; ++f) { 1770 PetscObject obj; 1771 PetscClassId id; 1772 1773 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1774 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1775 if (id == PETSCFV_CLASSID) useFVM = PETSC_TRUE; 1776 } 1777 /* Get local solution with boundary values */ 1778 ierr = DMGetLocalVector(dm, &locX);CHKERRQ(ierr); 1779 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 1780 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);CHKERRQ(ierr); 1781 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);CHKERRQ(ierr); 1782 /* Read DS information */ 1783 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1784 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 1785 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 1786 ierr = ISCreateStride(PETSC_COMM_SELF,numCells,cStart,1,&cellIS);CHKERRQ(ierr); 1787 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 1788 /* Read Auxiliary DS information */ 1789 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1790 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);CHKERRQ(ierr); 1791 if (dmAux) { 1792 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1793 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 1794 ierr = DMGetLocalSection(dmAux, §ionAux);CHKERRQ(ierr); 1795 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1796 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 1797 } 1798 /* Allocate data arrays */ 1799 ierr = PetscCalloc1(numCells*totDim, &u);CHKERRQ(ierr); 1800 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1801 /* Read out geometry */ 1802 ierr = DMGetCoordinateField(dm,&coordField);CHKERRQ(ierr); 1803 ierr = DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);CHKERRQ(ierr); 1804 if (maxDegree <= 1) { 1805 ierr = DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);CHKERRQ(ierr); 1806 if (affineQuad) { 1807 ierr = DMFieldCreateFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&cgeomFEM);CHKERRQ(ierr); 1808 } 1809 } 1810 if (useFVM) { 1811 PetscFV fv = NULL; 1812 Vec grad; 1813 PetscInt fStart, fEnd; 1814 PetscBool compGrad; 1815 1816 for (f = 0; f < Nf; ++f) { 1817 PetscObject obj; 1818 PetscClassId id; 1819 1820 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1821 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1822 if (id == PETSCFV_CLASSID) {fv = (PetscFV) obj; break;} 1823 } 1824 ierr = PetscFVGetComputeGradients(fv, &compGrad);CHKERRQ(ierr); 1825 ierr = PetscFVSetComputeGradients(fv, PETSC_TRUE);CHKERRQ(ierr); 1826 ierr = DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);CHKERRQ(ierr); 1827 ierr = DMPlexComputeGradientFVM(dm, fv, faceGeometryFVM, cellGeometryFVM, &dmGrad);CHKERRQ(ierr); 1828 ierr = PetscFVSetComputeGradients(fv, compGrad);CHKERRQ(ierr); 1829 ierr = VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1830 /* Reconstruct and limit cell gradients */ 1831 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1832 ierr = DMGetGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1833 ierr = DMPlexReconstructGradients_Internal(dm, fv, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);CHKERRQ(ierr); 1834 /* Communicate gradient values */ 1835 ierr = DMGetLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1836 ierr = DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1837 ierr = DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1838 ierr = DMRestoreGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1839 /* Handle non-essential (e.g. outflow) boundary values */ 1840 ierr = DMPlexInsertBoundaryValues(dm, PETSC_FALSE, locX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);CHKERRQ(ierr); 1841 ierr = VecGetArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1842 } 1843 /* Read out data from inputs */ 1844 for (c = cStart; c < cEnd; ++c) { 1845 PetscScalar *x = NULL; 1846 PetscInt i; 1847 1848 ierr = DMPlexVecGetClosure(dm, section, locX, c, NULL, &x);CHKERRQ(ierr); 1849 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 1850 ierr = DMPlexVecRestoreClosure(dm, section, locX, c, NULL, &x);CHKERRQ(ierr); 1851 if (dmAux) { 1852 ierr = DMPlexVecGetClosure(dmAux, sectionAux, locA, c, NULL, &x);CHKERRQ(ierr); 1853 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 1854 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, locA, c, NULL, &x);CHKERRQ(ierr); 1855 } 1856 } 1857 /* Do integration for each field */ 1858 for (f = 0; f < Nf; ++f) { 1859 PetscObject obj; 1860 PetscClassId id; 1861 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 1862 1863 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1864 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1865 if (id == PETSCFE_CLASSID) { 1866 PetscFE fe = (PetscFE) obj; 1867 PetscQuadrature q; 1868 PetscFEGeom *chunkGeom = NULL; 1869 PetscInt Nq, Nb; 1870 1871 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1872 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1873 ierr = PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 1874 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1875 blockSize = Nb*Nq; 1876 batchSize = numBlocks * blockSize; 1877 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1878 numChunks = numCells / (numBatches*batchSize); 1879 Ne = numChunks*numBatches*batchSize; 1880 Nr = numCells % (numBatches*batchSize); 1881 offset = numCells - Nr; 1882 if (!affineQuad) { 1883 ierr = DMFieldCreateFEGeom(coordField,cellIS,q,PETSC_FALSE,&cgeomFEM);CHKERRQ(ierr); 1884 } 1885 ierr = PetscFEGeomGetChunk(cgeomFEM,0,offset,&chunkGeom);CHKERRQ(ierr); 1886 ierr = PetscFEIntegrate(prob, f, Ne, chunkGeom, u, probAux, a, cintegral);CHKERRQ(ierr); 1887 ierr = PetscFEGeomGetChunk(cgeomFEM,offset,numCells,&chunkGeom);CHKERRQ(ierr); 1888 ierr = PetscFEIntegrate(prob, f, Nr, chunkGeom, &u[offset*totDim], probAux, &a[offset*totDimAux], &cintegral[offset*Nf]);CHKERRQ(ierr); 1889 ierr = PetscFEGeomRestoreChunk(cgeomFEM,offset,numCells,&chunkGeom);CHKERRQ(ierr); 1890 if (!affineQuad) { 1891 ierr = PetscFEGeomDestroy(&cgeomFEM);CHKERRQ(ierr); 1892 } 1893 } else if (id == PETSCFV_CLASSID) { 1894 PetscInt foff; 1895 PetscPointFunc obj_func; 1896 PetscScalar lint; 1897 1898 ierr = PetscDSGetObjective(prob, f, &obj_func);CHKERRQ(ierr); 1899 ierr = PetscDSGetFieldOffset(prob, f, &foff);CHKERRQ(ierr); 1900 if (obj_func) { 1901 for (c = 0; c < numCells; ++c) { 1902 PetscScalar *u_x; 1903 1904 ierr = DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);CHKERRQ(ierr); 1905 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); 1906 cintegral[c*Nf+f] += PetscRealPart(lint)*cgeomFVM[c].volume; 1907 } 1908 } 1909 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f); 1910 } 1911 /* Cleanup data arrays */ 1912 if (useFVM) { 1913 ierr = VecRestoreArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1914 ierr = VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1915 ierr = DMRestoreLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1916 ierr = VecDestroy(&faceGeometryFVM);CHKERRQ(ierr); 1917 ierr = VecDestroy(&cellGeometryFVM);CHKERRQ(ierr); 1918 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr); 1919 } 1920 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1921 ierr = PetscFree(u);CHKERRQ(ierr); 1922 /* Cleanup */ 1923 if (affineQuad) { 1924 ierr = PetscFEGeomDestroy(&cgeomFEM);CHKERRQ(ierr); 1925 } 1926 ierr = PetscQuadratureDestroy(&affineQuad);CHKERRQ(ierr); 1927 ierr = ISDestroy(&cellIS);CHKERRQ(ierr); 1928 ierr = DMRestoreLocalVector(dm, &locX);CHKERRQ(ierr); 1929 PetscFunctionReturn(0); 1930 } 1931 1932 /*@ 1933 DMPlexComputeIntegralFEM - Form the integral over the domain from the global input X using pointwise functions specified by the user 1934 1935 Input Parameters: 1936 + dm - The mesh 1937 . X - Global input vector 1938 - user - The user context 1939 1940 Output Parameter: 1941 . integral - Integral for each field 1942 1943 Level: developer 1944 1945 .seealso: DMPlexComputeResidualFEM() 1946 @*/ 1947 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscScalar *integral, void *user) 1948 { 1949 DM_Plex *mesh = (DM_Plex *) dm->data; 1950 PetscScalar *cintegral, *lintegral; 1951 PetscInt Nf, f, cellHeight, cStart, cEnd, cell; 1952 PetscErrorCode ierr; 1953 1954 PetscFunctionBegin; 1955 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1956 PetscValidHeaderSpecific(X, VEC_CLASSID, 2); 1957 PetscValidPointer(integral, 3); 1958 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1959 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 1960 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr); 1961 ierr = DMPlexGetSimplexOrBoxCells(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr); 1962 /* TODO Introduce a loop over large chunks (right now this is a single chunk) */ 1963 ierr = PetscCalloc2(Nf, &lintegral, (cEnd-cStart)*Nf, &cintegral);CHKERRQ(ierr); 1964 ierr = DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);CHKERRQ(ierr); 1965 /* Sum up values */ 1966 for (cell = cStart; cell < cEnd; ++cell) { 1967 const PetscInt c = cell - cStart; 1968 1969 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);CHKERRQ(ierr);} 1970 for (f = 0; f < Nf; ++f) lintegral[f] += cintegral[c*Nf+f]; 1971 } 1972 ierr = MPIU_Allreduce(lintegral, integral, Nf, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr); 1973 if (mesh->printFEM) { 1974 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Integral:");CHKERRQ(ierr); 1975 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", (double) PetscRealPart(integral[f]));CHKERRQ(ierr);} 1976 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 1977 } 1978 ierr = PetscFree2(lintegral, cintegral);CHKERRQ(ierr); 1979 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1980 PetscFunctionReturn(0); 1981 } 1982 1983 /*@ 1984 DMPlexComputeCellwiseIntegralFEM - Form the vector of cellwise integrals F from the global input X using pointwise functions specified by the user 1985 1986 Input Parameters: 1987 + dm - The mesh 1988 . X - Global input vector 1989 - user - The user context 1990 1991 Output Parameter: 1992 . integral - Cellwise integrals for each field 1993 1994 Level: developer 1995 1996 .seealso: DMPlexComputeResidualFEM() 1997 @*/ 1998 PetscErrorCode DMPlexComputeCellwiseIntegralFEM(DM dm, Vec X, Vec F, void *user) 1999 { 2000 DM_Plex *mesh = (DM_Plex *) dm->data; 2001 DM dmF; 2002 PetscSection sectionF; 2003 PetscScalar *cintegral, *af; 2004 PetscInt Nf, f, cellHeight, cStart, cEnd, cell; 2005 PetscErrorCode ierr; 2006 2007 PetscFunctionBegin; 2008 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2009 PetscValidHeaderSpecific(X, VEC_CLASSID, 2); 2010 PetscValidHeaderSpecific(F, VEC_CLASSID, 3); 2011 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 2012 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 2013 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr); 2014 ierr = DMPlexGetSimplexOrBoxCells(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr); 2015 /* TODO Introduce a loop over large chunks (right now this is a single chunk) */ 2016 ierr = PetscCalloc1((cEnd-cStart)*Nf, &cintegral);CHKERRQ(ierr); 2017 ierr = DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);CHKERRQ(ierr); 2018 /* Put values in F*/ 2019 ierr = VecGetDM(F, &dmF);CHKERRQ(ierr); 2020 ierr = DMGetLocalSection(dmF, §ionF);CHKERRQ(ierr); 2021 ierr = VecGetArray(F, &af);CHKERRQ(ierr); 2022 for (cell = cStart; cell < cEnd; ++cell) { 2023 const PetscInt c = cell - cStart; 2024 PetscInt dof, off; 2025 2026 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);CHKERRQ(ierr);} 2027 ierr = PetscSectionGetDof(sectionF, cell, &dof);CHKERRQ(ierr); 2028 ierr = PetscSectionGetOffset(sectionF, cell, &off);CHKERRQ(ierr); 2029 if (dof != Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cell dofs %D != %D", dof, Nf); 2030 for (f = 0; f < Nf; ++f) af[off+f] = cintegral[c*Nf+f]; 2031 } 2032 ierr = VecRestoreArray(F, &af);CHKERRQ(ierr); 2033 ierr = PetscFree(cintegral);CHKERRQ(ierr); 2034 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 2035 PetscFunctionReturn(0); 2036 } 2037 2038 static PetscErrorCode DMPlexComputeBdIntegral_Internal(DM dm, Vec locX, IS pointIS, 2039 void (*func)(PetscInt, PetscInt, PetscInt, 2040 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 2041 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 2042 PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 2043 PetscScalar *fintegral, void *user) 2044 { 2045 DM plex = NULL, plexA = NULL; 2046 DMEnclosureType encAux; 2047 PetscDS prob, probAux = NULL; 2048 PetscSection section, sectionAux = NULL; 2049 Vec locA = NULL; 2050 DMField coordField; 2051 PetscInt Nf, totDim, *uOff, *uOff_x; 2052 PetscInt NfAux = 0, totDimAux = 0, *aOff = NULL; 2053 PetscScalar *u, *a = NULL; 2054 const PetscScalar *constants; 2055 PetscInt numConstants, f; 2056 PetscErrorCode ierr; 2057 2058 PetscFunctionBegin; 2059 ierr = DMGetCoordinateField(dm, &coordField);CHKERRQ(ierr); 2060 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 2061 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 2062 ierr = DMGetLocalSection(dm, §ion);CHKERRQ(ierr); 2063 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 2064 /* Determine which discretizations we have */ 2065 for (f = 0; f < Nf; ++f) { 2066 PetscObject obj; 2067 PetscClassId id; 2068 2069 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 2070 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2071 if (id == PETSCFV_CLASSID) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Not supported for FVM (field %D)", f); 2072 } 2073 /* Read DS information */ 2074 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 2075 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 2076 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 2077 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 2078 /* Read Auxiliary DS information */ 2079 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);CHKERRQ(ierr); 2080 if (locA) { 2081 DM dmAux; 2082 2083 ierr = VecGetDM(locA, &dmAux);CHKERRQ(ierr); 2084 ierr = DMGetEnclosureRelation(dmAux, dm, &encAux);CHKERRQ(ierr); 2085 ierr = DMConvert(dmAux, DMPLEX, &plexA);CHKERRQ(ierr); 2086 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 2087 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 2088 ierr = DMGetLocalSection(dmAux, §ionAux);CHKERRQ(ierr); 2089 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 2090 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 2091 } 2092 /* Integrate over points */ 2093 { 2094 PetscFEGeom *fgeom, *chunkGeom = NULL; 2095 PetscInt maxDegree; 2096 PetscQuadrature qGeom = NULL; 2097 const PetscInt *points; 2098 PetscInt numFaces, face, Nq, field; 2099 PetscInt numChunks, chunkSize, chunk, Nr, offset; 2100 2101 ierr = ISGetLocalSize(pointIS, &numFaces);CHKERRQ(ierr); 2102 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 2103 ierr = PetscCalloc2(numFaces*totDim, &u, locA ? numFaces*totDimAux : 0, &a);CHKERRQ(ierr); 2104 ierr = DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree);CHKERRQ(ierr); 2105 for (field = 0; field < Nf; ++field) { 2106 PetscFE fe; 2107 2108 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr); 2109 if (maxDegree <= 1) {ierr = DMFieldCreateDefaultQuadrature(coordField, pointIS, &qGeom);CHKERRQ(ierr);} 2110 if (!qGeom) { 2111 ierr = PetscFEGetFaceQuadrature(fe, &qGeom);CHKERRQ(ierr); 2112 ierr = PetscObjectReference((PetscObject) qGeom);CHKERRQ(ierr); 2113 } 2114 ierr = PetscQuadratureGetData(qGeom, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 2115 ierr = DMPlexGetFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom);CHKERRQ(ierr); 2116 for (face = 0; face < numFaces; ++face) { 2117 const PetscInt point = points[face], *support, *cone; 2118 PetscScalar *x = NULL; 2119 PetscInt i, coneSize, faceLoc; 2120 2121 ierr = DMPlexGetSupport(dm, point, &support);CHKERRQ(ierr); 2122 ierr = DMPlexGetConeSize(dm, support[0], &coneSize);CHKERRQ(ierr); 2123 ierr = DMPlexGetCone(dm, support[0], &cone);CHKERRQ(ierr); 2124 for (faceLoc = 0; faceLoc < coneSize; ++faceLoc) if (cone[faceLoc] == point) break; 2125 if (faceLoc == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %D in cone of support[0] %D", face, support[0]); 2126 fgeom->face[face][0] = faceLoc; 2127 ierr = DMPlexVecGetClosure(plex, section, locX, support[0], NULL, &x);CHKERRQ(ierr); 2128 for (i = 0; i < totDim; ++i) u[face*totDim+i] = x[i]; 2129 ierr = DMPlexVecRestoreClosure(plex, section, locX, support[0], NULL, &x);CHKERRQ(ierr); 2130 if (locA) { 2131 PetscInt subp; 2132 ierr = DMGetEnclosurePoint(plexA, dm, encAux, support[0], &subp);CHKERRQ(ierr); 2133 ierr = DMPlexVecGetClosure(plexA, sectionAux, locA, subp, NULL, &x);CHKERRQ(ierr); 2134 for (i = 0; i < totDimAux; ++i) a[f*totDimAux+i] = x[i]; 2135 ierr = DMPlexVecRestoreClosure(plexA, sectionAux, locA, subp, NULL, &x);CHKERRQ(ierr); 2136 } 2137 } 2138 /* Get blocking */ 2139 { 2140 PetscQuadrature q; 2141 PetscInt numBatches, batchSize, numBlocks, blockSize; 2142 PetscInt Nq, Nb; 2143 2144 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 2145 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 2146 ierr = PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 2147 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 2148 blockSize = Nb*Nq; 2149 batchSize = numBlocks * blockSize; 2150 chunkSize = numBatches*batchSize; 2151 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 2152 numChunks = numFaces / chunkSize; 2153 Nr = numFaces % chunkSize; 2154 offset = numFaces - Nr; 2155 } 2156 /* Do integration for each field */ 2157 for (chunk = 0; chunk < numChunks; ++chunk) { 2158 ierr = PetscFEGeomGetChunk(fgeom, chunk*chunkSize, (chunk+1)*chunkSize, &chunkGeom);CHKERRQ(ierr); 2159 ierr = PetscFEIntegrateBd(prob, field, func, chunkSize, chunkGeom, u, probAux, a, fintegral);CHKERRQ(ierr); 2160 ierr = PetscFEGeomRestoreChunk(fgeom, 0, offset, &chunkGeom);CHKERRQ(ierr); 2161 } 2162 ierr = PetscFEGeomGetChunk(fgeom, offset, numFaces, &chunkGeom);CHKERRQ(ierr); 2163 ierr = PetscFEIntegrateBd(prob, field, func, Nr, chunkGeom, &u[offset*totDim], probAux, a ? &a[offset*totDimAux] : NULL, &fintegral[offset*Nf]);CHKERRQ(ierr); 2164 ierr = PetscFEGeomRestoreChunk(fgeom, offset, numFaces, &chunkGeom);CHKERRQ(ierr); 2165 /* Cleanup data arrays */ 2166 ierr = DMPlexRestoreFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom);CHKERRQ(ierr); 2167 ierr = PetscQuadratureDestroy(&qGeom);CHKERRQ(ierr); 2168 ierr = PetscFree2(u, a);CHKERRQ(ierr); 2169 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 2170 } 2171 } 2172 if (plex) {ierr = DMDestroy(&plex);CHKERRQ(ierr);} 2173 if (plexA) {ierr = DMDestroy(&plexA);CHKERRQ(ierr);} 2174 PetscFunctionReturn(0); 2175 } 2176 2177 /*@ 2178 DMPlexComputeBdIntegral - Form the integral over the specified boundary from the global input X using pointwise functions specified by the user 2179 2180 Input Parameters: 2181 + dm - The mesh 2182 . X - Global input vector 2183 . label - The boundary DMLabel 2184 . numVals - The number of label values to use, or PETSC_DETERMINE for all values 2185 . vals - The label values to use, or PETSC_NULL for all values 2186 . func = The function to integrate along the boundary 2187 - user - The user context 2188 2189 Output Parameter: 2190 . integral - Integral for each field 2191 2192 Level: developer 2193 2194 .seealso: DMPlexComputeIntegralFEM(), DMPlexComputeBdResidualFEM() 2195 @*/ 2196 PetscErrorCode DMPlexComputeBdIntegral(DM dm, Vec X, DMLabel label, PetscInt numVals, const PetscInt vals[], 2197 void (*func)(PetscInt, PetscInt, PetscInt, 2198 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 2199 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 2200 PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), 2201 PetscScalar *integral, void *user) 2202 { 2203 Vec locX; 2204 PetscSection section; 2205 DMLabel depthLabel; 2206 IS facetIS; 2207 PetscInt dim, Nf, f, v; 2208 PetscErrorCode ierr; 2209 2210 PetscFunctionBegin; 2211 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2212 PetscValidHeaderSpecific(X, VEC_CLASSID, 2); 2213 PetscValidPointer(label, 3); 2214 if (vals) PetscValidPointer(vals, 5); 2215 PetscValidPointer(integral, 6); 2216 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 2217 ierr = DMPlexGetDepthLabel(dm, &depthLabel);CHKERRQ(ierr); 2218 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2219 ierr = DMLabelGetStratumIS(depthLabel, dim-1, &facetIS);CHKERRQ(ierr); 2220 ierr = DMGetLocalSection(dm, §ion);CHKERRQ(ierr); 2221 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 2222 /* Get local solution with boundary values */ 2223 ierr = DMGetLocalVector(dm, &locX);CHKERRQ(ierr); 2224 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 2225 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);CHKERRQ(ierr); 2226 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);CHKERRQ(ierr); 2227 /* Loop over label values */ 2228 ierr = PetscArrayzero(integral, Nf);CHKERRQ(ierr); 2229 for (v = 0; v < numVals; ++v) { 2230 IS pointIS; 2231 PetscInt numFaces, face; 2232 PetscScalar *fintegral; 2233 2234 ierr = DMLabelGetStratumIS(label, vals[v], &pointIS);CHKERRQ(ierr); 2235 if (!pointIS) continue; /* No points with that id on this process */ 2236 { 2237 IS isectIS; 2238 2239 /* TODO: Special cases of ISIntersect where it is quick to check a priori if one is a superset of the other */ 2240 ierr = ISIntersect_Caching_Internal(facetIS, pointIS, &isectIS);CHKERRQ(ierr); 2241 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 2242 pointIS = isectIS; 2243 } 2244 ierr = ISGetLocalSize(pointIS, &numFaces);CHKERRQ(ierr); 2245 ierr = PetscCalloc1(numFaces*Nf, &fintegral);CHKERRQ(ierr); 2246 ierr = DMPlexComputeBdIntegral_Internal(dm, locX, pointIS, func, fintegral, user);CHKERRQ(ierr); 2247 /* Sum point contributions into integral */ 2248 for (f = 0; f < Nf; ++f) for (face = 0; face < numFaces; ++face) integral[f] += fintegral[face*Nf+f]; 2249 ierr = PetscFree(fintegral);CHKERRQ(ierr); 2250 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 2251 } 2252 ierr = DMRestoreLocalVector(dm, &locX);CHKERRQ(ierr); 2253 ierr = ISDestroy(&facetIS);CHKERRQ(ierr); 2254 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 2255 PetscFunctionReturn(0); 2256 } 2257 2258 /*@ 2259 DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to a uniformly refined DM. 2260 2261 Input Parameters: 2262 + dmc - The coarse mesh 2263 . dmf - The fine mesh 2264 . isRefined - Flag indicating regular refinement, rather than the same topology 2265 - user - The user context 2266 2267 Output Parameter: 2268 . In - The interpolation matrix 2269 2270 Level: developer 2271 2272 .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 2273 @*/ 2274 PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, PetscBool isRefined, Mat In, void *user) 2275 { 2276 DM_Plex *mesh = (DM_Plex *) dmc->data; 2277 const char *name = "Interpolator"; 2278 PetscDS cds, rds; 2279 PetscFE *feRef; 2280 PetscFV *fvRef; 2281 PetscSection fsection, fglobalSection; 2282 PetscSection csection, cglobalSection; 2283 PetscScalar *elemMat; 2284 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, c; 2285 PetscInt cTotDim, rTotDim = 0; 2286 PetscErrorCode ierr; 2287 2288 PetscFunctionBegin; 2289 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2290 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 2291 ierr = DMGetLocalSection(dmf, &fsection);CHKERRQ(ierr); 2292 ierr = DMGetGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 2293 ierr = DMGetLocalSection(dmc, &csection);CHKERRQ(ierr); 2294 ierr = DMGetGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 2295 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 2296 ierr = DMPlexGetSimplexOrBoxCells(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 2297 ierr = DMGetDS(dmc, &cds);CHKERRQ(ierr); 2298 ierr = DMGetDS(dmf, &rds);CHKERRQ(ierr); 2299 ierr = PetscCalloc2(Nf, &feRef, Nf, &fvRef);CHKERRQ(ierr); 2300 for (f = 0; f < Nf; ++f) { 2301 PetscObject obj; 2302 PetscClassId id; 2303 PetscInt rNb = 0, Nc = 0; 2304 2305 ierr = PetscDSGetDiscretization(rds, f, &obj);CHKERRQ(ierr); 2306 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2307 if (id == PETSCFE_CLASSID) { 2308 PetscFE fe = (PetscFE) obj; 2309 2310 if (isRefined) { 2311 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 2312 } else { 2313 ierr = PetscObjectReference((PetscObject) fe);CHKERRQ(ierr); 2314 feRef[f] = fe; 2315 } 2316 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 2317 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 2318 } else if (id == PETSCFV_CLASSID) { 2319 PetscFV fv = (PetscFV) obj; 2320 PetscDualSpace Q; 2321 2322 if (isRefined) { 2323 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 2324 } else { 2325 ierr = PetscObjectReference((PetscObject) fv);CHKERRQ(ierr); 2326 fvRef[f] = fv; 2327 } 2328 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 2329 ierr = PetscDualSpaceGetDimension(Q, &rNb);CHKERRQ(ierr); 2330 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 2331 } 2332 rTotDim += rNb; 2333 } 2334 ierr = PetscDSGetTotalDimension(cds, &cTotDim);CHKERRQ(ierr); 2335 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 2336 ierr = PetscArrayzero(elemMat, rTotDim*cTotDim);CHKERRQ(ierr); 2337 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 2338 PetscDualSpace Qref; 2339 PetscQuadrature f; 2340 const PetscReal *qpoints, *qweights; 2341 PetscReal *points; 2342 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 2343 2344 /* Compose points from all dual basis functionals */ 2345 if (feRef[fieldI]) { 2346 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 2347 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 2348 } else { 2349 ierr = PetscFVGetDualSpace(fvRef[fieldI], &Qref);CHKERRQ(ierr); 2350 ierr = PetscFVGetNumComponents(fvRef[fieldI], &Nc);CHKERRQ(ierr); 2351 } 2352 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 2353 for (i = 0; i < fpdim; ++i) { 2354 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 2355 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 2356 npoints += Np; 2357 } 2358 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 2359 for (i = 0, k = 0; i < fpdim; ++i) { 2360 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 2361 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 2362 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 2363 } 2364 2365 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 2366 PetscObject obj; 2367 PetscClassId id; 2368 PetscInt NcJ = 0, cpdim = 0, j, qNc; 2369 2370 ierr = PetscDSGetDiscretization(cds, fieldJ, &obj);CHKERRQ(ierr); 2371 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2372 if (id == PETSCFE_CLASSID) { 2373 PetscFE fe = (PetscFE) obj; 2374 PetscTabulation T = NULL; 2375 2376 /* Evaluate basis at points */ 2377 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 2378 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 2379 /* For now, fields only interpolate themselves */ 2380 if (fieldI == fieldJ) { 2381 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); 2382 ierr = PetscFECreateTabulation(fe, 1, npoints, points, 0, &T);CHKERRQ(ierr); 2383 for (i = 0, k = 0; i < fpdim; ++i) { 2384 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 2385 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);CHKERRQ(ierr); 2386 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); 2387 for (p = 0; p < Np; ++p, ++k) { 2388 for (j = 0; j < cpdim; ++j) { 2389 /* 2390 cTotDim: Total columns in element interpolation matrix, sum of number of dual basis functionals in each field 2391 offsetI, offsetJ: Offsets into the larger element interpolation matrix for different fields 2392 fpdim, i, cpdim, j: Dofs for fine and coarse grids, correspond to dual space basis functionals 2393 qNC, Nc, Ncj, c: Number of components in this field 2394 Np, p: Number of quad points in the fine grid functional i 2395 k: i*Np + p, overall point number for the interpolation 2396 */ 2397 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]; 2398 } 2399 } 2400 } 2401 ierr = PetscTabulationDestroy(&T);CHKERRQ(ierr);CHKERRQ(ierr); 2402 } 2403 } else if (id == PETSCFV_CLASSID) { 2404 PetscFV fv = (PetscFV) obj; 2405 2406 /* Evaluate constant function at points */ 2407 ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr); 2408 cpdim = 1; 2409 /* For now, fields only interpolate themselves */ 2410 if (fieldI == fieldJ) { 2411 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); 2412 for (i = 0, k = 0; i < fpdim; ++i) { 2413 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 2414 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);CHKERRQ(ierr); 2415 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); 2416 for (p = 0; p < Np; ++p, ++k) { 2417 for (j = 0; j < cpdim; ++j) { 2418 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += 1.0*qweights[p*qNc+c]; 2419 } 2420 } 2421 } 2422 } 2423 } 2424 offsetJ += cpdim; 2425 } 2426 offsetI += fpdim; 2427 ierr = PetscFree(points);CHKERRQ(ierr); 2428 } 2429 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 2430 /* Preallocate matrix */ 2431 { 2432 Mat preallocator; 2433 PetscScalar *vals; 2434 PetscInt *cellCIndices, *cellFIndices; 2435 PetscInt locRows, locCols, cell; 2436 2437 ierr = MatGetLocalSize(In, &locRows, &locCols);CHKERRQ(ierr); 2438 ierr = MatCreate(PetscObjectComm((PetscObject) In), &preallocator);CHKERRQ(ierr); 2439 ierr = MatSetType(preallocator, MATPREALLOCATOR);CHKERRQ(ierr); 2440 ierr = MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 2441 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 2442 ierr = PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 2443 for (cell = cStart; cell < cEnd; ++cell) { 2444 if (isRefined) { 2445 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 2446 ierr = MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);CHKERRQ(ierr); 2447 } else { 2448 ierr = DMPlexMatSetClosureGeneral(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, preallocator, cell, vals, INSERT_VALUES);CHKERRQ(ierr); 2449 } 2450 } 2451 ierr = PetscFree3(vals,cellCIndices,cellFIndices);CHKERRQ(ierr); 2452 ierr = MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2453 ierr = MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2454 ierr = MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);CHKERRQ(ierr); 2455 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 2456 } 2457 /* Fill matrix */ 2458 ierr = MatZeroEntries(In);CHKERRQ(ierr); 2459 for (c = cStart; c < cEnd; ++c) { 2460 if (isRefined) { 2461 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 2462 } else { 2463 ierr = DMPlexMatSetClosureGeneral(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 2464 } 2465 } 2466 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 2467 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 2468 ierr = PetscFree(elemMat);CHKERRQ(ierr); 2469 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2470 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2471 if (mesh->printFEM) { 2472 ierr = PetscPrintf(PetscObjectComm((PetscObject)In), "%s:\n", name);CHKERRQ(ierr); 2473 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 2474 ierr = MatView(In, NULL);CHKERRQ(ierr); 2475 } 2476 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2477 PetscFunctionReturn(0); 2478 } 2479 2480 PetscErrorCode DMPlexComputeMassMatrixNested(DM dmc, DM dmf, Mat mass, void *user) 2481 { 2482 SETERRQ(PetscObjectComm((PetscObject) dmc), PETSC_ERR_SUP, "Laziness"); 2483 } 2484 2485 /*@ 2486 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 2487 2488 Input Parameters: 2489 + dmf - The fine mesh 2490 . dmc - The coarse mesh 2491 - user - The user context 2492 2493 Output Parameter: 2494 . In - The interpolation matrix 2495 2496 Level: developer 2497 2498 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 2499 @*/ 2500 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 2501 { 2502 DM_Plex *mesh = (DM_Plex *) dmf->data; 2503 const char *name = "Interpolator"; 2504 PetscDS prob; 2505 PetscSection fsection, csection, globalFSection, globalCSection; 2506 PetscHSetIJ ht; 2507 PetscLayout rLayout; 2508 PetscInt *dnz, *onz; 2509 PetscInt locRows, rStart, rEnd; 2510 PetscReal *x, *v0, *J, *invJ, detJ; 2511 PetscReal *v0c, *Jc, *invJc, detJc; 2512 PetscScalar *elemMat; 2513 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 2514 PetscErrorCode ierr; 2515 2516 PetscFunctionBegin; 2517 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2518 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 2519 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 2520 ierr = PetscDSGetWorkspace(prob, &x, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 2521 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 2522 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 2523 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 2524 ierr = DMGetLocalSection(dmf, &fsection);CHKERRQ(ierr); 2525 ierr = DMGetGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 2526 ierr = DMGetLocalSection(dmc, &csection);CHKERRQ(ierr); 2527 ierr = DMGetGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 2528 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 2529 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 2530 ierr = PetscMalloc1(totDim, &elemMat);CHKERRQ(ierr); 2531 2532 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 2533 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 2534 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 2535 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 2536 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 2537 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 2538 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 2539 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 2540 ierr = PetscHSetIJCreate(&ht);CHKERRQ(ierr); 2541 for (field = 0; field < Nf; ++field) { 2542 PetscObject obj; 2543 PetscClassId id; 2544 PetscDualSpace Q = NULL; 2545 PetscQuadrature f; 2546 const PetscReal *qpoints; 2547 PetscInt Nc, Np, fpdim, i, d; 2548 2549 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 2550 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2551 if (id == PETSCFE_CLASSID) { 2552 PetscFE fe = (PetscFE) obj; 2553 2554 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 2555 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 2556 } else if (id == PETSCFV_CLASSID) { 2557 PetscFV fv = (PetscFV) obj; 2558 2559 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 2560 Nc = 1; 2561 } 2562 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 2563 /* For each fine grid cell */ 2564 for (cell = cStart; cell < cEnd; ++cell) { 2565 PetscInt *findices, *cindices; 2566 PetscInt numFIndices, numCIndices; 2567 2568 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);CHKERRQ(ierr); 2569 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 2570 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %D != %D dual basis vecs", numFIndices, fpdim); 2571 for (i = 0; i < fpdim; ++i) { 2572 Vec pointVec; 2573 PetscScalar *pV; 2574 PetscSF coarseCellSF = NULL; 2575 const PetscSFNode *coarseCells; 2576 PetscInt numCoarseCells, q, c; 2577 2578 /* Get points from the dual basis functional quadrature */ 2579 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 2580 ierr = PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 2581 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 2582 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 2583 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2584 for (q = 0; q < Np; ++q) { 2585 const PetscReal xi0[3] = {-1., -1., -1.}; 2586 2587 /* Transform point to real space */ 2588 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x); 2589 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 2590 } 2591 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2592 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 2593 /* OPT: Pack all quad points from fine cell */ 2594 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 2595 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 2596 /* Update preallocation info */ 2597 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 2598 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 2599 { 2600 PetscHashIJKey key; 2601 PetscBool missing; 2602 2603 key.i = findices[i]; 2604 if (key.i >= 0) { 2605 /* Get indices for coarse elements */ 2606 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 2607 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);CHKERRQ(ierr); 2608 for (c = 0; c < numCIndices; ++c) { 2609 key.j = cindices[c]; 2610 if (key.j < 0) continue; 2611 ierr = PetscHSetIJQueryAdd(ht, key, &missing);CHKERRQ(ierr); 2612 if (missing) { 2613 if ((key.j >= rStart) && (key.j < rEnd)) ++dnz[key.i-rStart]; 2614 else ++onz[key.i-rStart]; 2615 } 2616 } 2617 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);CHKERRQ(ierr); 2618 } 2619 } 2620 } 2621 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 2622 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 2623 } 2624 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);CHKERRQ(ierr); 2625 } 2626 } 2627 ierr = PetscHSetIJDestroy(&ht);CHKERRQ(ierr); 2628 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 2629 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2630 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 2631 for (field = 0; field < Nf; ++field) { 2632 PetscObject obj; 2633 PetscClassId id; 2634 PetscDualSpace Q = NULL; 2635 PetscTabulation T = NULL; 2636 PetscQuadrature f; 2637 const PetscReal *qpoints, *qweights; 2638 PetscInt Nc, qNc, Np, fpdim, i, d; 2639 2640 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 2641 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2642 if (id == PETSCFE_CLASSID) { 2643 PetscFE fe = (PetscFE) obj; 2644 2645 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 2646 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 2647 ierr = PetscFECreateTabulation(fe, 1, 1, x, 0, &T);CHKERRQ(ierr); 2648 } else if (id == PETSCFV_CLASSID) { 2649 PetscFV fv = (PetscFV) obj; 2650 2651 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 2652 Nc = 1; 2653 } else SETERRQ1(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONG,"Unknown discretization type for field %D",field); 2654 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 2655 /* For each fine grid cell */ 2656 for (cell = cStart; cell < cEnd; ++cell) { 2657 PetscInt *findices, *cindices; 2658 PetscInt numFIndices, numCIndices; 2659 2660 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);CHKERRQ(ierr); 2661 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 2662 if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %D != %D dual basis vecs", numFIndices, fpdim); 2663 for (i = 0; i < fpdim; ++i) { 2664 Vec pointVec; 2665 PetscScalar *pV; 2666 PetscSF coarseCellSF = NULL; 2667 const PetscSFNode *coarseCells; 2668 PetscInt numCoarseCells, cpdim, q, c, j; 2669 2670 /* Get points from the dual basis functional quadrature */ 2671 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 2672 ierr = PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);CHKERRQ(ierr); 2673 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); 2674 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 2675 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 2676 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2677 for (q = 0; q < Np; ++q) { 2678 const PetscReal xi0[3] = {-1., -1., -1.}; 2679 2680 /* Transform point to real space */ 2681 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x); 2682 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 2683 } 2684 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2685 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 2686 /* OPT: Read this out from preallocation information */ 2687 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 2688 /* Update preallocation info */ 2689 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 2690 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 2691 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2692 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 2693 PetscReal pVReal[3]; 2694 const PetscReal xi0[3] = {-1., -1., -1.}; 2695 2696 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);CHKERRQ(ierr); 2697 /* Transform points from real space to coarse reference space */ 2698 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 2699 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 2700 CoordinatesRealToRef(dim, dim, xi0, v0c, invJc, pVReal, x); 2701 2702 if (id == PETSCFE_CLASSID) { 2703 PetscFE fe = (PetscFE) obj; 2704 2705 /* Evaluate coarse basis on contained point */ 2706 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 2707 ierr = PetscFEComputeTabulation(fe, 1, x, 0, T);CHKERRQ(ierr); 2708 ierr = PetscArrayzero(elemMat, cpdim);CHKERRQ(ierr); 2709 /* Get elemMat entries by multiplying by weight */ 2710 for (j = 0; j < cpdim; ++j) { 2711 for (c = 0; c < Nc; ++c) elemMat[j] += T->T[0][j*Nc + c]*qweights[ccell*qNc + c]; 2712 } 2713 } else { 2714 cpdim = 1; 2715 for (j = 0; j < cpdim; ++j) { 2716 for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c]; 2717 } 2718 } 2719 /* Update interpolator */ 2720 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 2721 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 2722 ierr = MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 2723 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);CHKERRQ(ierr); 2724 } 2725 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2726 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 2727 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 2728 } 2729 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);CHKERRQ(ierr); 2730 } 2731 if (id == PETSCFE_CLASSID) {ierr = PetscTabulationDestroy(&T);CHKERRQ(ierr);} 2732 } 2733 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 2734 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 2735 ierr = PetscFree(elemMat);CHKERRQ(ierr); 2736 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2737 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2738 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 2739 PetscFunctionReturn(0); 2740 } 2741 2742 /*@ 2743 DMPlexComputeMassMatrixGeneral - Form the local portion of the mass matrix M from the coarse DM to a non-nested fine DM. 2744 2745 Input Parameters: 2746 + dmf - The fine mesh 2747 . dmc - The coarse mesh 2748 - user - The user context 2749 2750 Output Parameter: 2751 . mass - The mass matrix 2752 2753 Level: developer 2754 2755 .seealso: DMPlexComputeMassMatrixNested(), DMPlexComputeInterpolatorNested(), DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 2756 @*/ 2757 PetscErrorCode DMPlexComputeMassMatrixGeneral(DM dmc, DM dmf, Mat mass, void *user) 2758 { 2759 DM_Plex *mesh = (DM_Plex *) dmf->data; 2760 const char *name = "Mass Matrix"; 2761 PetscDS prob; 2762 PetscSection fsection, csection, globalFSection, globalCSection; 2763 PetscHSetIJ ht; 2764 PetscLayout rLayout; 2765 PetscInt *dnz, *onz; 2766 PetscInt locRows, rStart, rEnd; 2767 PetscReal *x, *v0, *J, *invJ, detJ; 2768 PetscReal *v0c, *Jc, *invJc, detJc; 2769 PetscScalar *elemMat; 2770 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 2771 PetscErrorCode ierr; 2772 2773 PetscFunctionBegin; 2774 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 2775 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 2776 ierr = PetscDSGetWorkspace(prob, &x, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 2777 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 2778 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 2779 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 2780 ierr = DMGetLocalSection(dmf, &fsection);CHKERRQ(ierr); 2781 ierr = DMGetGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 2782 ierr = DMGetLocalSection(dmc, &csection);CHKERRQ(ierr); 2783 ierr = DMGetGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 2784 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 2785 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 2786 ierr = PetscMalloc1(totDim, &elemMat);CHKERRQ(ierr); 2787 2788 ierr = MatGetLocalSize(mass, &locRows, NULL);CHKERRQ(ierr); 2789 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) mass), &rLayout);CHKERRQ(ierr); 2790 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 2791 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 2792 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 2793 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 2794 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 2795 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 2796 ierr = PetscHSetIJCreate(&ht);CHKERRQ(ierr); 2797 for (field = 0; field < Nf; ++field) { 2798 PetscObject obj; 2799 PetscClassId id; 2800 PetscQuadrature quad; 2801 const PetscReal *qpoints; 2802 PetscInt Nq, Nc, i, d; 2803 2804 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 2805 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2806 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetQuadrature((PetscFE) obj, &quad);CHKERRQ(ierr);} 2807 else {ierr = PetscFVGetQuadrature((PetscFV) obj, &quad);CHKERRQ(ierr);} 2808 ierr = PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, NULL);CHKERRQ(ierr); 2809 /* For each fine grid cell */ 2810 for (cell = cStart; cell < cEnd; ++cell) { 2811 Vec pointVec; 2812 PetscScalar *pV; 2813 PetscSF coarseCellSF = NULL; 2814 const PetscSFNode *coarseCells; 2815 PetscInt numCoarseCells, q, c; 2816 PetscInt *findices, *cindices; 2817 PetscInt numFIndices, numCIndices; 2818 2819 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);CHKERRQ(ierr); 2820 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 2821 /* Get points from the quadrature */ 2822 ierr = VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);CHKERRQ(ierr); 2823 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 2824 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2825 for (q = 0; q < Nq; ++q) { 2826 const PetscReal xi0[3] = {-1., -1., -1.}; 2827 2828 /* Transform point to real space */ 2829 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x); 2830 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 2831 } 2832 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2833 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 2834 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 2835 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 2836 /* Update preallocation info */ 2837 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 2838 if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 2839 { 2840 PetscHashIJKey key; 2841 PetscBool missing; 2842 2843 for (i = 0; i < numFIndices; ++i) { 2844 key.i = findices[i]; 2845 if (key.i >= 0) { 2846 /* Get indices for coarse elements */ 2847 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 2848 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);CHKERRQ(ierr); 2849 for (c = 0; c < numCIndices; ++c) { 2850 key.j = cindices[c]; 2851 if (key.j < 0) continue; 2852 ierr = PetscHSetIJQueryAdd(ht, key, &missing);CHKERRQ(ierr); 2853 if (missing) { 2854 if ((key.j >= rStart) && (key.j < rEnd)) ++dnz[key.i-rStart]; 2855 else ++onz[key.i-rStart]; 2856 } 2857 } 2858 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);CHKERRQ(ierr); 2859 } 2860 } 2861 } 2862 } 2863 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 2864 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 2865 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);CHKERRQ(ierr); 2866 } 2867 } 2868 ierr = PetscHSetIJDestroy(&ht);CHKERRQ(ierr); 2869 ierr = MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 2870 ierr = MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2871 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 2872 for (field = 0; field < Nf; ++field) { 2873 PetscObject obj; 2874 PetscClassId id; 2875 PetscTabulation T, Tfine; 2876 PetscQuadrature quad; 2877 const PetscReal *qpoints, *qweights; 2878 PetscInt Nq, Nc, i, d; 2879 2880 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 2881 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 2882 if (id == PETSCFE_CLASSID) { 2883 ierr = PetscFEGetQuadrature((PetscFE) obj, &quad);CHKERRQ(ierr); 2884 ierr = PetscFEGetCellTabulation((PetscFE) obj, &Tfine);CHKERRQ(ierr); 2885 ierr = PetscFECreateTabulation((PetscFE) obj, 1, 1, x, 0, &T);CHKERRQ(ierr); 2886 } else { 2887 ierr = PetscFVGetQuadrature((PetscFV) obj, &quad);CHKERRQ(ierr); 2888 } 2889 ierr = PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, &qweights);CHKERRQ(ierr); 2890 /* For each fine grid cell */ 2891 for (cell = cStart; cell < cEnd; ++cell) { 2892 Vec pointVec; 2893 PetscScalar *pV; 2894 PetscSF coarseCellSF = NULL; 2895 const PetscSFNode *coarseCells; 2896 PetscInt numCoarseCells, cpdim, q, c, j; 2897 PetscInt *findices, *cindices; 2898 PetscInt numFIndices, numCIndices; 2899 2900 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);CHKERRQ(ierr); 2901 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 2902 /* Get points from the quadrature */ 2903 ierr = VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);CHKERRQ(ierr); 2904 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 2905 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2906 for (q = 0; q < Nq; ++q) { 2907 const PetscReal xi0[3] = {-1., -1., -1.}; 2908 2909 /* Transform point to real space */ 2910 CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x); 2911 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 2912 } 2913 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2914 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 2915 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 2916 /* Update matrix */ 2917 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 2918 if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 2919 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 2920 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 2921 PetscReal pVReal[3]; 2922 const PetscReal xi0[3] = {-1., -1., -1.}; 2923 2924 2925 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);CHKERRQ(ierr); 2926 /* Transform points from real space to coarse reference space */ 2927 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 2928 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 2929 CoordinatesRealToRef(dim, dim, xi0, v0c, invJc, pVReal, x); 2930 2931 if (id == PETSCFE_CLASSID) { 2932 PetscFE fe = (PetscFE) obj; 2933 2934 /* Evaluate coarse basis on contained point */ 2935 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 2936 ierr = PetscFEComputeTabulation(fe, 1, x, 0, T);CHKERRQ(ierr); 2937 /* Get elemMat entries by multiplying by weight */ 2938 for (i = 0; i < numFIndices; ++i) { 2939 ierr = PetscArrayzero(elemMat, cpdim);CHKERRQ(ierr); 2940 for (j = 0; j < cpdim; ++j) { 2941 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; 2942 } 2943 /* Update interpolator */ 2944 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 2945 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 2946 ierr = MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);CHKERRQ(ierr); 2947 } 2948 } else { 2949 cpdim = 1; 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] += 1.0*1.0*qweights[ccell*Nc + c]*detJ; 2954 } 2955 /* Update interpolator */ 2956 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);} 2957 ierr = PetscPrintf(PETSC_COMM_SELF, "Nq: %D %D Nf: %D %D Nc: %D %D\n", ccell, Nq, i, numFIndices, j, numCIndices);CHKERRQ(ierr); 2958 if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim); 2959 ierr = MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);CHKERRQ(ierr); 2960 } 2961 } 2962 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);CHKERRQ(ierr); 2963 } 2964 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 2965 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 2966 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 2967 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);CHKERRQ(ierr); 2968 } 2969 if (id == PETSCFE_CLASSID) {ierr = PetscTabulationDestroy(&T);CHKERRQ(ierr);} 2970 } 2971 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 2972 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 2973 ierr = PetscFree(elemMat);CHKERRQ(ierr); 2974 ierr = MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2975 ierr = MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2976 PetscFunctionReturn(0); 2977 } 2978 2979 /*@ 2980 DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns 2981 2982 Input Parameters: 2983 + dmc - The coarse mesh 2984 - dmf - The fine mesh 2985 - user - The user context 2986 2987 Output Parameter: 2988 . sc - The mapping 2989 2990 Level: developer 2991 2992 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 2993 @*/ 2994 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 2995 { 2996 PetscDS prob; 2997 PetscFE *feRef; 2998 PetscFV *fvRef; 2999 Vec fv, cv; 3000 IS fis, cis; 3001 PetscSection fsection, fglobalSection, csection, cglobalSection; 3002 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 3003 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, c, dim, d, startC, endC, offsetC, offsetF, m; 3004 PetscBool *needAvg; 3005 PetscErrorCode ierr; 3006 3007 PetscFunctionBegin; 3008 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 3009 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 3010 ierr = DMGetLocalSection(dmf, &fsection);CHKERRQ(ierr); 3011 ierr = DMGetGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 3012 ierr = DMGetLocalSection(dmc, &csection);CHKERRQ(ierr); 3013 ierr = DMGetGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 3014 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 3015 ierr = DMPlexGetSimplexOrBoxCells(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 3016 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 3017 ierr = PetscCalloc3(Nf,&feRef,Nf,&fvRef,Nf,&needAvg);CHKERRQ(ierr); 3018 for (f = 0; f < Nf; ++f) { 3019 PetscObject obj; 3020 PetscClassId id; 3021 PetscInt fNb = 0, Nc = 0; 3022 3023 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 3024 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 3025 if (id == PETSCFE_CLASSID) { 3026 PetscFE fe = (PetscFE) obj; 3027 PetscSpace sp; 3028 PetscInt maxDegree; 3029 3030 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 3031 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 3032 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 3033 ierr = PetscFEGetBasisSpace(fe, &sp);CHKERRQ(ierr); 3034 ierr = PetscSpaceGetDegree(sp, NULL, &maxDegree);CHKERRQ(ierr); 3035 if (!maxDegree) needAvg[f] = PETSC_TRUE; 3036 } else if (id == PETSCFV_CLASSID) { 3037 PetscFV fv = (PetscFV) obj; 3038 PetscDualSpace Q; 3039 3040 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 3041 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 3042 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 3043 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 3044 needAvg[f] = PETSC_TRUE; 3045 } 3046 fTotDim += fNb; 3047 } 3048 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 3049 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 3050 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 3051 PetscFE feC; 3052 PetscFV fvC; 3053 PetscDualSpace QF, QC; 3054 PetscInt order = -1, NcF, NcC, fpdim, cpdim; 3055 3056 if (feRef[field]) { 3057 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 3058 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 3059 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 3060 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 3061 ierr = PetscDualSpaceGetOrder(QF, &order);CHKERRQ(ierr); 3062 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 3063 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 3064 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 3065 } else { 3066 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 3067 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 3068 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 3069 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 3070 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 3071 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 3072 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 3073 } 3074 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); 3075 for (c = 0; c < cpdim; ++c) { 3076 PetscQuadrature cfunc; 3077 const PetscReal *cqpoints, *cqweights; 3078 PetscInt NqcC, NpC; 3079 PetscBool found = PETSC_FALSE; 3080 3081 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 3082 ierr = PetscQuadratureGetData(cfunc, NULL, &NqcC, &NpC, &cqpoints, &cqweights);CHKERRQ(ierr); 3083 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); 3084 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 3085 for (f = 0; f < fpdim; ++f) { 3086 PetscQuadrature ffunc; 3087 const PetscReal *fqpoints, *fqweights; 3088 PetscReal sum = 0.0; 3089 PetscInt NqcF, NpF; 3090 3091 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 3092 ierr = PetscQuadratureGetData(ffunc, NULL, &NqcF, &NpF, &fqpoints, &fqweights);CHKERRQ(ierr); 3093 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); 3094 if (NpC != NpF) continue; 3095 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 3096 if (sum > 1.0e-9) continue; 3097 for (d = 0; d < NcC; ++d) sum += PetscAbsReal(cqweights[d]*fqweights[d]); 3098 if (sum < 1.0e-9) continue; 3099 cmap[offsetC+c] = offsetF+f; 3100 found = PETSC_TRUE; 3101 break; 3102 } 3103 if (!found) { 3104 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 3105 if (fvRef[field] || (feRef[field] && order == 0)) { 3106 cmap[offsetC+c] = offsetF+0; 3107 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 3108 } 3109 } 3110 offsetC += cpdim; 3111 offsetF += fpdim; 3112 } 3113 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 3114 ierr = PetscFree3(feRef,fvRef,needAvg);CHKERRQ(ierr); 3115 3116 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 3117 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 3118 ierr = VecGetOwnershipRange(cv, &startC, &endC);CHKERRQ(ierr); 3119 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 3120 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 3121 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 3122 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 3123 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 3124 for (c = cStart; c < cEnd; ++c) { 3125 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 3126 for (d = 0; d < cTotDim; ++d) { 3127 if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue; 3128 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]]); 3129 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 3130 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 3131 } 3132 } 3133 ierr = PetscFree(cmap);CHKERRQ(ierr); 3134 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 3135 3136 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 3137 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 3138 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 3139 ierr = ISDestroy(&cis);CHKERRQ(ierr); 3140 ierr = ISDestroy(&fis);CHKERRQ(ierr); 3141 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 3142 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 3143 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 3144 PetscFunctionReturn(0); 3145 } 3146 3147 /*@C 3148 DMPlexGetCellFields - Retrieve the field values values for a chunk of cells 3149 3150 Input Parameters: 3151 + dm - The DM 3152 . cellIS - The cells to include 3153 . locX - A local vector with the solution fields 3154 . locX_t - A local vector with solution field time derivatives, or NULL 3155 - locA - A local vector with auxiliary fields, or NULL 3156 3157 Output Parameters: 3158 + u - The field coefficients 3159 . u_t - The fields derivative coefficients 3160 - a - The auxiliary field coefficients 3161 3162 Level: developer 3163 3164 .seealso: DMPlexGetFaceFields() 3165 @*/ 3166 PetscErrorCode DMPlexGetCellFields(DM dm, IS cellIS, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a) 3167 { 3168 DM plex, plexA = NULL; 3169 DMEnclosureType encAux; 3170 PetscSection section, sectionAux; 3171 PetscDS prob; 3172 const PetscInt *cells; 3173 PetscInt cStart, cEnd, numCells, totDim, totDimAux, c; 3174 PetscErrorCode ierr; 3175 3176 PetscFunctionBegin; 3177 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3178 PetscValidHeaderSpecific(locX, VEC_CLASSID, 4); 3179 if (locX_t) {PetscValidHeaderSpecific(locX_t, VEC_CLASSID, 5);} 3180 if (locA) {PetscValidHeaderSpecific(locA, VEC_CLASSID, 6);} 3181 PetscValidPointer(u, 7); 3182 PetscValidPointer(u_t, 8); 3183 PetscValidPointer(a, 9); 3184 ierr = DMPlexConvertPlex(dm, &plex, PETSC_FALSE);CHKERRQ(ierr); 3185 ierr = ISGetPointRange(cellIS, &cStart, &cEnd, &cells);CHKERRQ(ierr); 3186 ierr = DMGetLocalSection(dm, §ion);CHKERRQ(ierr); 3187 ierr = DMGetCellDS(dm, cStart, &prob);CHKERRQ(ierr); 3188 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 3189 if (locA) { 3190 DM dmAux; 3191 PetscDS probAux; 3192 3193 ierr = VecGetDM(locA, &dmAux);CHKERRQ(ierr); 3194 ierr = DMGetEnclosureRelation(dmAux, dm, &encAux);CHKERRQ(ierr); 3195 ierr = DMPlexConvertPlex(dmAux, &plexA, PETSC_FALSE);CHKERRQ(ierr); 3196 ierr = DMGetLocalSection(dmAux, §ionAux);CHKERRQ(ierr); 3197 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 3198 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 3199 } 3200 numCells = cEnd - cStart; 3201 ierr = DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, u);CHKERRQ(ierr); 3202 if (locX_t) {ierr = DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, u_t);CHKERRQ(ierr);} else {*u_t = NULL;} 3203 if (locA) {ierr = DMGetWorkArray(dm, numCells*totDimAux, MPIU_SCALAR, a);CHKERRQ(ierr);} else {*a = NULL;} 3204 for (c = cStart; c < cEnd; ++c) { 3205 const PetscInt cell = cells ? cells[c] : c; 3206 const PetscInt cind = c - cStart; 3207 PetscScalar *x = NULL, *x_t = NULL, *ul = *u, *ul_t = *u_t, *al = *a; 3208 PetscInt i; 3209 3210 ierr = DMPlexVecGetClosure(plex, section, locX, cell, NULL, &x);CHKERRQ(ierr); 3211 for (i = 0; i < totDim; ++i) ul[cind*totDim+i] = x[i]; 3212 ierr = DMPlexVecRestoreClosure(plex, section, locX, cell, NULL, &x);CHKERRQ(ierr); 3213 if (locX_t) { 3214 ierr = DMPlexVecGetClosure(plex, section, locX_t, cell, NULL, &x_t);CHKERRQ(ierr); 3215 for (i = 0; i < totDim; ++i) ul_t[cind*totDim+i] = x_t[i]; 3216 ierr = DMPlexVecRestoreClosure(plex, section, locX_t, cell, NULL, &x_t);CHKERRQ(ierr); 3217 } 3218 if (locA) { 3219 PetscInt subcell; 3220 ierr = DMGetEnclosurePoint(plexA, dm, encAux, cell, &subcell);CHKERRQ(ierr); 3221 ierr = DMPlexVecGetClosure(plexA, sectionAux, locA, subcell, NULL, &x);CHKERRQ(ierr); 3222 for (i = 0; i < totDimAux; ++i) al[cind*totDimAux+i] = x[i]; 3223 ierr = DMPlexVecRestoreClosure(plexA, sectionAux, locA, subcell, NULL, &x);CHKERRQ(ierr); 3224 } 3225 } 3226 ierr = DMDestroy(&plex);CHKERRQ(ierr); 3227 if (locA) {ierr = DMDestroy(&plexA);CHKERRQ(ierr);} 3228 ierr = ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);CHKERRQ(ierr); 3229 PetscFunctionReturn(0); 3230 } 3231 3232 /*@C 3233 DMPlexRestoreCellFields - Restore the field values values for a chunk of cells 3234 3235 Input Parameters: 3236 + dm - The DM 3237 . cellIS - The cells to include 3238 . locX - A local vector with the solution fields 3239 . locX_t - A local vector with solution field time derivatives, or NULL 3240 - locA - A local vector with auxiliary fields, or NULL 3241 3242 Output Parameters: 3243 + u - The field coefficients 3244 . u_t - The fields derivative coefficients 3245 - a - The auxiliary field coefficients 3246 3247 Level: developer 3248 3249 .seealso: DMPlexGetFaceFields() 3250 @*/ 3251 PetscErrorCode DMPlexRestoreCellFields(DM dm, IS cellIS, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a) 3252 { 3253 PetscErrorCode ierr; 3254 3255 PetscFunctionBegin; 3256 ierr = DMRestoreWorkArray(dm, 0, MPIU_SCALAR, u);CHKERRQ(ierr); 3257 if (locX_t) {ierr = DMRestoreWorkArray(dm, 0, MPIU_SCALAR, u_t);CHKERRQ(ierr);} 3258 if (locA) {ierr = DMRestoreWorkArray(dm, 0, MPIU_SCALAR, a);CHKERRQ(ierr);} 3259 PetscFunctionReturn(0); 3260 } 3261 3262 /*@C 3263 DMPlexGetFaceFields - Retrieve the field values values for a chunk of faces 3264 3265 Input Parameters: 3266 + dm - The DM 3267 . fStart - The first face to include 3268 . fEnd - The first face to exclude 3269 . locX - A local vector with the solution fields 3270 . locX_t - A local vector with solution field time derivatives, or NULL 3271 . faceGeometry - A local vector with face geometry 3272 . cellGeometry - A local vector with cell geometry 3273 - locaGrad - A local vector with field gradients, or NULL 3274 3275 Output Parameters: 3276 + Nface - The number of faces with field values 3277 . uL - The field values at the left side of the face 3278 - uR - The field values at the right side of the face 3279 3280 Level: developer 3281 3282 .seealso: DMPlexGetCellFields() 3283 @*/ 3284 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) 3285 { 3286 DM dmFace, dmCell, dmGrad = NULL; 3287 PetscSection section; 3288 PetscDS prob; 3289 DMLabel ghostLabel; 3290 const PetscScalar *facegeom, *cellgeom, *x, *lgrad; 3291 PetscBool *isFE; 3292 PetscInt dim, Nf, f, Nc, numFaces = fEnd - fStart, iface, face; 3293 PetscErrorCode ierr; 3294 3295 PetscFunctionBegin; 3296 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3297 PetscValidHeaderSpecific(locX, VEC_CLASSID, 4); 3298 if (locX_t) {PetscValidHeaderSpecific(locX_t, VEC_CLASSID, 5);} 3299 PetscValidHeaderSpecific(faceGeometry, VEC_CLASSID, 6); 3300 PetscValidHeaderSpecific(cellGeometry, VEC_CLASSID, 7); 3301 if (locGrad) {PetscValidHeaderSpecific(locGrad, VEC_CLASSID, 8);} 3302 PetscValidPointer(uL, 9); 3303 PetscValidPointer(uR, 10); 3304 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3305 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 3306 ierr = DMGetLocalSection(dm, §ion);CHKERRQ(ierr); 3307 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 3308 ierr = PetscDSGetTotalComponents(prob, &Nc);CHKERRQ(ierr); 3309 ierr = PetscMalloc1(Nf, &isFE);CHKERRQ(ierr); 3310 for (f = 0; f < Nf; ++f) { 3311 PetscObject obj; 3312 PetscClassId id; 3313 3314 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 3315 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 3316 if (id == PETSCFE_CLASSID) {isFE[f] = PETSC_TRUE;} 3317 else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;} 3318 else {isFE[f] = PETSC_FALSE;} 3319 } 3320 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 3321 ierr = VecGetArrayRead(locX, &x);CHKERRQ(ierr); 3322 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 3323 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 3324 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 3325 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 3326 if (locGrad) { 3327 ierr = VecGetDM(locGrad, &dmGrad);CHKERRQ(ierr); 3328 ierr = VecGetArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 3329 } 3330 ierr = DMGetWorkArray(dm, numFaces*Nc, MPIU_SCALAR, uL);CHKERRQ(ierr); 3331 ierr = DMGetWorkArray(dm, numFaces*Nc, MPIU_SCALAR, uR);CHKERRQ(ierr); 3332 /* Right now just eat the extra work for FE (could make a cell loop) */ 3333 for (face = fStart, iface = 0; face < fEnd; ++face) { 3334 const PetscInt *cells; 3335 PetscFVFaceGeom *fg; 3336 PetscFVCellGeom *cgL, *cgR; 3337 PetscScalar *xL, *xR, *gL, *gR; 3338 PetscScalar *uLl = *uL, *uRl = *uR; 3339 PetscInt ghost, nsupp, nchild; 3340 3341 ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr); 3342 ierr = DMPlexGetSupportSize(dm, face, &nsupp);CHKERRQ(ierr); 3343 ierr = DMPlexGetTreeChildren(dm, face, &nchild, NULL);CHKERRQ(ierr); 3344 if (ghost >= 0 || nsupp > 2 || nchild > 0) continue; 3345 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 3346 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 3347 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);CHKERRQ(ierr); 3348 ierr = DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);CHKERRQ(ierr); 3349 for (f = 0; f < Nf; ++f) { 3350 PetscInt off; 3351 3352 ierr = PetscDSGetComponentOffset(prob, f, &off);CHKERRQ(ierr); 3353 if (isFE[f]) { 3354 const PetscInt *cone; 3355 PetscInt comp, coneSizeL, coneSizeR, faceLocL, faceLocR, ldof, rdof, d; 3356 3357 xL = xR = NULL; 3358 ierr = PetscSectionGetFieldComponents(section, f, &comp);CHKERRQ(ierr); 3359 ierr = DMPlexVecGetClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);CHKERRQ(ierr); 3360 ierr = DMPlexVecGetClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);CHKERRQ(ierr); 3361 ierr = DMPlexGetCone(dm, cells[0], &cone);CHKERRQ(ierr); 3362 ierr = DMPlexGetConeSize(dm, cells[0], &coneSizeL);CHKERRQ(ierr); 3363 for (faceLocL = 0; faceLocL < coneSizeL; ++faceLocL) if (cone[faceLocL] == face) break; 3364 ierr = DMPlexGetCone(dm, cells[1], &cone);CHKERRQ(ierr); 3365 ierr = DMPlexGetConeSize(dm, cells[1], &coneSizeR);CHKERRQ(ierr); 3366 for (faceLocR = 0; faceLocR < coneSizeR; ++faceLocR) if (cone[faceLocR] == face) break; 3367 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]); 3368 /* Check that FEM field has values in the right cell (sometimes its an FV ghost cell) */ 3369 /* TODO: this is a hack that might not be right for nonconforming */ 3370 if (faceLocL < coneSizeL) { 3371 ierr = PetscFEEvaluateFaceFields_Internal(prob, f, faceLocL, xL, &uLl[iface*Nc+off]);CHKERRQ(ierr); 3372 if (rdof == ldof && faceLocR < coneSizeR) {ierr = PetscFEEvaluateFaceFields_Internal(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);CHKERRQ(ierr);} 3373 else {for(d = 0; d < comp; ++d) uRl[iface*Nc+off+d] = uLl[iface*Nc+off+d];} 3374 } 3375 else { 3376 ierr = PetscFEEvaluateFaceFields_Internal(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);CHKERRQ(ierr); 3377 ierr = PetscSectionGetFieldComponents(section, f, &comp);CHKERRQ(ierr); 3378 for(d = 0; d < comp; ++d) uLl[iface*Nc+off+d] = uRl[iface*Nc+off+d]; 3379 } 3380 ierr = DMPlexVecRestoreClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);CHKERRQ(ierr); 3381 ierr = DMPlexVecRestoreClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);CHKERRQ(ierr); 3382 } else { 3383 PetscFV fv; 3384 PetscInt numComp, c; 3385 3386 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fv);CHKERRQ(ierr); 3387 ierr = PetscFVGetNumComponents(fv, &numComp);CHKERRQ(ierr); 3388 ierr = DMPlexPointLocalFieldRead(dm, cells[0], f, x, &xL);CHKERRQ(ierr); 3389 ierr = DMPlexPointLocalFieldRead(dm, cells[1], f, x, &xR);CHKERRQ(ierr); 3390 if (dmGrad) { 3391 PetscReal dxL[3], dxR[3]; 3392 3393 ierr = DMPlexPointLocalRead(dmGrad, cells[0], lgrad, &gL);CHKERRQ(ierr); 3394 ierr = DMPlexPointLocalRead(dmGrad, cells[1], lgrad, &gR);CHKERRQ(ierr); 3395 DMPlex_WaxpyD_Internal(dim, -1, cgL->centroid, fg->centroid, dxL); 3396 DMPlex_WaxpyD_Internal(dim, -1, cgR->centroid, fg->centroid, dxR); 3397 for (c = 0; c < numComp; ++c) { 3398 uLl[iface*Nc+off+c] = xL[c] + DMPlex_DotD_Internal(dim, &gL[c*dim], dxL); 3399 uRl[iface*Nc+off+c] = xR[c] + DMPlex_DotD_Internal(dim, &gR[c*dim], dxR); 3400 } 3401 } else { 3402 for (c = 0; c < numComp; ++c) { 3403 uLl[iface*Nc+off+c] = xL[c]; 3404 uRl[iface*Nc+off+c] = xR[c]; 3405 } 3406 } 3407 } 3408 } 3409 ++iface; 3410 } 3411 *Nface = iface; 3412 ierr = VecRestoreArrayRead(locX, &x);CHKERRQ(ierr); 3413 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 3414 ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 3415 if (locGrad) { 3416 ierr = VecRestoreArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 3417 } 3418 ierr = PetscFree(isFE);CHKERRQ(ierr); 3419 PetscFunctionReturn(0); 3420 } 3421 3422 /*@C 3423 DMPlexRestoreFaceFields - Restore the field values values for a chunk of faces 3424 3425 Input Parameters: 3426 + dm - The DM 3427 . fStart - The first face to include 3428 . fEnd - The first face to exclude 3429 . locX - A local vector with the solution fields 3430 . locX_t - A local vector with solution field time derivatives, or NULL 3431 . faceGeometry - A local vector with face geometry 3432 . cellGeometry - A local vector with cell geometry 3433 - locaGrad - A local vector with field gradients, or NULL 3434 3435 Output Parameters: 3436 + Nface - The number of faces with field values 3437 . uL - The field values at the left side of the face 3438 - uR - The field values at the right side of the face 3439 3440 Level: developer 3441 3442 .seealso: DMPlexGetFaceFields() 3443 @*/ 3444 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) 3445 { 3446 PetscErrorCode ierr; 3447 3448 PetscFunctionBegin; 3449 ierr = DMRestoreWorkArray(dm, 0, MPIU_SCALAR, uL);CHKERRQ(ierr); 3450 ierr = DMRestoreWorkArray(dm, 0, MPIU_SCALAR, uR);CHKERRQ(ierr); 3451 PetscFunctionReturn(0); 3452 } 3453 3454 /*@C 3455 DMPlexGetFaceGeometry - Retrieve the geometric values for a chunk of faces 3456 3457 Input Parameters: 3458 + dm - The DM 3459 . fStart - The first face to include 3460 . fEnd - The first face to exclude 3461 . faceGeometry - A local vector with face geometry 3462 - cellGeometry - A local vector with cell geometry 3463 3464 Output Parameters: 3465 + Nface - The number of faces with field values 3466 . fgeom - The extract the face centroid and normal 3467 - vol - The cell volume 3468 3469 Level: developer 3470 3471 .seealso: DMPlexGetCellFields() 3472 @*/ 3473 PetscErrorCode DMPlexGetFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscInt *Nface, PetscFVFaceGeom **fgeom, PetscReal **vol) 3474 { 3475 DM dmFace, dmCell; 3476 DMLabel ghostLabel; 3477 const PetscScalar *facegeom, *cellgeom; 3478 PetscInt dim, numFaces = fEnd - fStart, iface, face; 3479 PetscErrorCode ierr; 3480 3481 PetscFunctionBegin; 3482 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3483 PetscValidHeaderSpecific(faceGeometry, VEC_CLASSID, 4); 3484 PetscValidHeaderSpecific(cellGeometry, VEC_CLASSID, 5); 3485 PetscValidPointer(fgeom, 6); 3486 PetscValidPointer(vol, 7); 3487 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 3488 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 3489 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 3490 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 3491 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 3492 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 3493 ierr = PetscMalloc1(numFaces, fgeom);CHKERRQ(ierr); 3494 ierr = DMGetWorkArray(dm, numFaces*2, MPIU_SCALAR, vol);CHKERRQ(ierr); 3495 for (face = fStart, iface = 0; face < fEnd; ++face) { 3496 const PetscInt *cells; 3497 PetscFVFaceGeom *fg; 3498 PetscFVCellGeom *cgL, *cgR; 3499 PetscFVFaceGeom *fgeoml = *fgeom; 3500 PetscReal *voll = *vol; 3501 PetscInt ghost, d, nchild, nsupp; 3502 3503 ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr); 3504 ierr = DMPlexGetSupportSize(dm, face, &nsupp);CHKERRQ(ierr); 3505 ierr = DMPlexGetTreeChildren(dm, face, &nchild, NULL);CHKERRQ(ierr); 3506 if (ghost >= 0 || nsupp > 2 || nchild > 0) continue; 3507 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 3508 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 3509 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);CHKERRQ(ierr); 3510 ierr = DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);CHKERRQ(ierr); 3511 for (d = 0; d < dim; ++d) { 3512 fgeoml[iface].centroid[d] = fg->centroid[d]; 3513 fgeoml[iface].normal[d] = fg->normal[d]; 3514 } 3515 voll[iface*2+0] = cgL->volume; 3516 voll[iface*2+1] = cgR->volume; 3517 ++iface; 3518 } 3519 *Nface = iface; 3520 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 3521 ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 3522 PetscFunctionReturn(0); 3523 } 3524 3525 /*@C 3526 DMPlexRestoreFaceGeometry - Restore the field values values for a chunk of faces 3527 3528 Input Parameters: 3529 + dm - The DM 3530 . fStart - The first face to include 3531 . fEnd - The first face to exclude 3532 . faceGeometry - A local vector with face geometry 3533 - cellGeometry - A local vector with cell geometry 3534 3535 Output Parameters: 3536 + Nface - The number of faces with field values 3537 . fgeom - The extract the face centroid and normal 3538 - vol - The cell volume 3539 3540 Level: developer 3541 3542 .seealso: DMPlexGetFaceFields() 3543 @*/ 3544 PetscErrorCode DMPlexRestoreFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscInt *Nface, PetscFVFaceGeom **fgeom, PetscReal **vol) 3545 { 3546 PetscErrorCode ierr; 3547 3548 PetscFunctionBegin; 3549 ierr = PetscFree(*fgeom);CHKERRQ(ierr); 3550 ierr = DMRestoreWorkArray(dm, 0, MPIU_REAL, vol);CHKERRQ(ierr); 3551 PetscFunctionReturn(0); 3552 } 3553 3554 PetscErrorCode DMSNESGetFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom) 3555 { 3556 char composeStr[33] = {0}; 3557 PetscObjectId id; 3558 PetscContainer container; 3559 PetscErrorCode ierr; 3560 3561 PetscFunctionBegin; 3562 ierr = PetscObjectGetId((PetscObject)quad,&id);CHKERRQ(ierr); 3563 ierr = PetscSNPrintf(composeStr, 32, "DMSNESGetFEGeom_%x\n", id);CHKERRQ(ierr); 3564 ierr = PetscObjectQuery((PetscObject) pointIS, composeStr, (PetscObject *) &container);CHKERRQ(ierr); 3565 if (container) { 3566 ierr = PetscContainerGetPointer(container, (void **) geom);CHKERRQ(ierr); 3567 } else { 3568 ierr = DMFieldCreateFEGeom(coordField, pointIS, quad, faceData, geom);CHKERRQ(ierr); 3569 ierr = PetscContainerCreate(PETSC_COMM_SELF,&container);CHKERRQ(ierr); 3570 ierr = PetscContainerSetPointer(container, (void *) *geom);CHKERRQ(ierr); 3571 ierr = PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom);CHKERRQ(ierr); 3572 ierr = PetscObjectCompose((PetscObject) pointIS, composeStr, (PetscObject) container);CHKERRQ(ierr); 3573 ierr = PetscContainerDestroy(&container);CHKERRQ(ierr); 3574 } 3575 PetscFunctionReturn(0); 3576 } 3577 3578 PetscErrorCode DMSNESRestoreFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom) 3579 { 3580 PetscFunctionBegin; 3581 *geom = NULL; 3582 PetscFunctionReturn(0); 3583 } 3584 3585 PetscErrorCode DMPlexComputeResidual_Patch_Internal(DM dm, PetscSection section, IS cellIS, PetscReal t, Vec locX, Vec locX_t, Vec locF, void *user) 3586 { 3587 DM_Plex *mesh = (DM_Plex *) dm->data; 3588 const char *name = "Residual"; 3589 DM dmAux = NULL; 3590 DMLabel ghostLabel = NULL; 3591 PetscDS prob = NULL; 3592 PetscDS probAux = NULL; 3593 PetscBool useFEM = PETSC_FALSE; 3594 PetscBool isImplicit = (locX_t || t == PETSC_MIN_REAL) ? PETSC_TRUE : PETSC_FALSE; 3595 DMField coordField = NULL; 3596 Vec locA; 3597 PetscScalar *u = NULL, *u_t, *a, *uL = NULL, *uR = NULL; 3598 IS chunkIS; 3599 const PetscInt *cells; 3600 PetscInt cStart, cEnd, numCells; 3601 PetscInt Nf, f, totDim, totDimAux, numChunks, cellChunkSize, chunk, fStart, fEnd; 3602 PetscInt maxDegree = PETSC_MAX_INT; 3603 PetscQuadrature affineQuad = NULL, *quads = NULL; 3604 PetscFEGeom *affineGeom = NULL, **geoms = NULL; 3605 PetscErrorCode ierr; 3606 3607 PetscFunctionBegin; 3608 ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr); 3609 /* FEM+FVM */ 3610 /* 1: Get sizes from dm and dmAux */ 3611 ierr = DMGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 3612 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 3613 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 3614 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 3615 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);CHKERRQ(ierr); 3616 if (locA) { 3617 ierr = VecGetDM(locA, &dmAux);CHKERRQ(ierr); 3618 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 3619 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 3620 } 3621 /* 2: Get geometric data */ 3622 for (f = 0; f < Nf; ++f) { 3623 PetscObject obj; 3624 PetscClassId id; 3625 PetscBool fimp; 3626 3627 ierr = PetscDSGetImplicit(prob, f, &fimp);CHKERRQ(ierr); 3628 if (isImplicit != fimp) continue; 3629 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 3630 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 3631 if (id == PETSCFE_CLASSID) {useFEM = PETSC_TRUE;} 3632 if (id == PETSCFV_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Use of FVM with PCPATCH not yet implemented"); 3633 } 3634 if (useFEM) { 3635 ierr = DMGetCoordinateField(dm, &coordField);CHKERRQ(ierr); 3636 ierr = DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);CHKERRQ(ierr); 3637 if (maxDegree <= 1) { 3638 ierr = DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);CHKERRQ(ierr); 3639 if (affineQuad) { 3640 ierr = DMSNESGetFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);CHKERRQ(ierr); 3641 } 3642 } else { 3643 ierr = PetscCalloc2(Nf,&quads,Nf,&geoms);CHKERRQ(ierr); 3644 for (f = 0; f < Nf; ++f) { 3645 PetscObject obj; 3646 PetscClassId id; 3647 PetscBool fimp; 3648 3649 ierr = PetscDSGetImplicit(prob, f, &fimp);CHKERRQ(ierr); 3650 if (isImplicit != fimp) continue; 3651 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 3652 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 3653 if (id == PETSCFE_CLASSID) { 3654 PetscFE fe = (PetscFE) obj; 3655 3656 ierr = PetscFEGetQuadrature(fe, &quads[f]);CHKERRQ(ierr); 3657 ierr = PetscObjectReference((PetscObject)quads[f]);CHKERRQ(ierr); 3658 ierr = DMSNESGetFEGeom(coordField,cellIS,quads[f],PETSC_FALSE,&geoms[f]);CHKERRQ(ierr); 3659 } 3660 } 3661 } 3662 } 3663 /* Loop over chunks */ 3664 ierr = ISGetPointRange(cellIS, &cStart, &cEnd, &cells);CHKERRQ(ierr); 3665 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 3666 if (useFEM) {ierr = ISCreate(PETSC_COMM_SELF, &chunkIS);CHKERRQ(ierr);} 3667 numCells = cEnd - cStart; 3668 numChunks = 1; 3669 cellChunkSize = numCells/numChunks; 3670 numChunks = PetscMin(1,numCells); 3671 for (chunk = 0; chunk < numChunks; ++chunk) { 3672 PetscScalar *elemVec, *fluxL = NULL, *fluxR = NULL; 3673 PetscReal *vol = NULL; 3674 PetscFVFaceGeom *fgeom = NULL; 3675 PetscInt cS = cStart+chunk*cellChunkSize, cE = PetscMin(cS+cellChunkSize, cEnd), numCells = cE - cS, c; 3676 PetscInt numFaces = 0; 3677 3678 /* Extract field coefficients */ 3679 if (useFEM) { 3680 ierr = ISGetPointSubrange(chunkIS, cS, cE, cells);CHKERRQ(ierr); 3681 ierr = DMPlexGetCellFields(dm, chunkIS, locX, locX_t, locA, &u, &u_t, &a);CHKERRQ(ierr); 3682 ierr = DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);CHKERRQ(ierr); 3683 ierr = PetscArrayzero(elemVec, numCells*totDim);CHKERRQ(ierr); 3684 } 3685 /* TODO We will interlace both our field coefficients (u, u_t, uL, uR, etc.) and our output (elemVec, fL, fR). I think this works */ 3686 /* Loop over fields */ 3687 for (f = 0; f < Nf; ++f) { 3688 PetscObject obj; 3689 PetscClassId id; 3690 PetscBool fimp; 3691 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 3692 3693 ierr = PetscDSGetImplicit(prob, f, &fimp);CHKERRQ(ierr); 3694 if (isImplicit != fimp) continue; 3695 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 3696 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 3697 if (id == PETSCFE_CLASSID) { 3698 PetscFE fe = (PetscFE) obj; 3699 PetscFEGeom *geom = affineGeom ? affineGeom : geoms[f]; 3700 PetscFEGeom *chunkGeom = NULL; 3701 PetscQuadrature quad = affineQuad ? affineQuad : quads[f]; 3702 PetscInt Nq, Nb; 3703 3704 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 3705 ierr = PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 3706 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 3707 blockSize = Nb; 3708 batchSize = numBlocks * blockSize; 3709 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 3710 numChunks = numCells / (numBatches*batchSize); 3711 Ne = numChunks*numBatches*batchSize; 3712 Nr = numCells % (numBatches*batchSize); 3713 offset = numCells - Nr; 3714 /* Integrate FE residual to get elemVec (need fields at quadrature points) */ 3715 /* 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) */ 3716 ierr = PetscFEGeomGetChunk(geom,0,offset,&chunkGeom);CHKERRQ(ierr); 3717 ierr = PetscFEIntegrateResidual(prob, f, Ne, chunkGeom, u, u_t, probAux, a, t, elemVec);CHKERRQ(ierr); 3718 ierr = PetscFEGeomGetChunk(geom,offset,numCells,&chunkGeom);CHKERRQ(ierr); 3719 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); 3720 ierr = PetscFEGeomRestoreChunk(geom,offset,numCells,&chunkGeom);CHKERRQ(ierr); 3721 } else if (id == PETSCFV_CLASSID) { 3722 PetscFV fv = (PetscFV) obj; 3723 3724 Ne = numFaces; 3725 /* Riemann solve over faces (need fields at face centroids) */ 3726 /* We need to evaluate FE fields at those coordinates */ 3727 ierr = PetscFVIntegrateRHSFunction(fv, prob, f, Ne, fgeom, vol, uL, uR, fluxL, fluxR);CHKERRQ(ierr); 3728 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f); 3729 } 3730 /* Loop over domain */ 3731 if (useFEM) { 3732 /* Add elemVec to locX */ 3733 for (c = cS; c < cE; ++c) { 3734 const PetscInt cell = cells ? cells[c] : c; 3735 const PetscInt cind = c - cStart; 3736 3737 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(cell, name, totDim, &elemVec[cind*totDim]);CHKERRQ(ierr);} 3738 if (ghostLabel) { 3739 PetscInt ghostVal; 3740 3741 ierr = DMLabelGetValue(ghostLabel,cell,&ghostVal);CHKERRQ(ierr); 3742 if (ghostVal > 0) continue; 3743 } 3744 ierr = DMPlexVecSetClosure(dm, section, locF, cell, &elemVec[cind*totDim], ADD_ALL_VALUES);CHKERRQ(ierr); 3745 } 3746 } 3747 /* Handle time derivative */ 3748 if (locX_t) { 3749 PetscScalar *x_t, *fa; 3750 3751 ierr = VecGetArray(locF, &fa);CHKERRQ(ierr); 3752 ierr = VecGetArray(locX_t, &x_t);CHKERRQ(ierr); 3753 for (f = 0; f < Nf; ++f) { 3754 PetscFV fv; 3755 PetscObject obj; 3756 PetscClassId id; 3757 PetscInt pdim, d; 3758 3759 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 3760 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 3761 if (id != PETSCFV_CLASSID) continue; 3762 fv = (PetscFV) obj; 3763 ierr = PetscFVGetNumComponents(fv, &pdim);CHKERRQ(ierr); 3764 for (c = cS; c < cE; ++c) { 3765 const PetscInt cell = cells ? cells[c] : c; 3766 PetscScalar *u_t, *r; 3767 3768 if (ghostLabel) { 3769 PetscInt ghostVal; 3770 3771 ierr = DMLabelGetValue(ghostLabel, cell, &ghostVal);CHKERRQ(ierr); 3772 if (ghostVal > 0) continue; 3773 } 3774 ierr = DMPlexPointLocalFieldRead(dm, cell, f, x_t, &u_t);CHKERRQ(ierr); 3775 ierr = DMPlexPointLocalFieldRef(dm, cell, f, fa, &r);CHKERRQ(ierr); 3776 for (d = 0; d < pdim; ++d) r[d] += u_t[d]; 3777 } 3778 } 3779 ierr = VecRestoreArray(locX_t, &x_t);CHKERRQ(ierr); 3780 ierr = VecRestoreArray(locF, &fa);CHKERRQ(ierr); 3781 } 3782 if (useFEM) { 3783 ierr = DMPlexRestoreCellFields(dm, chunkIS, locX, locX_t, locA, &u, &u_t, &a);CHKERRQ(ierr); 3784 ierr = DMRestoreWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);CHKERRQ(ierr); 3785 } 3786 } 3787 if (useFEM) {ierr = ISDestroy(&chunkIS);CHKERRQ(ierr);} 3788 ierr = ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);CHKERRQ(ierr); 3789 /* TODO Could include boundary residual here (see DMPlexComputeResidual_Internal) */ 3790 if (useFEM) { 3791 if (maxDegree <= 1) { 3792 ierr = DMSNESRestoreFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);CHKERRQ(ierr); 3793 ierr = PetscQuadratureDestroy(&affineQuad);CHKERRQ(ierr); 3794 } else { 3795 for (f = 0; f < Nf; ++f) { 3796 ierr = DMSNESRestoreFEGeom(coordField,cellIS,quads[f],PETSC_FALSE,&geoms[f]);CHKERRQ(ierr); 3797 ierr = PetscQuadratureDestroy(&quads[f]);CHKERRQ(ierr); 3798 } 3799 ierr = PetscFree2(quads,geoms);CHKERRQ(ierr); 3800 } 3801 } 3802 ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr); 3803 PetscFunctionReturn(0); 3804 } 3805 3806 /* 3807 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 3808 3809 X - The local solution vector 3810 X_t - The local solution time derviative vector, or NULL 3811 */ 3812 PetscErrorCode DMPlexComputeJacobian_Patch_Internal(DM dm, PetscSection section, PetscSection globalSection, IS cellIS, 3813 PetscReal t, PetscReal X_tShift, Vec X, Vec X_t, Mat Jac, Mat JacP, void *ctx) 3814 { 3815 DM_Plex *mesh = (DM_Plex *) dm->data; 3816 const char *name = "Jacobian", *nameP = "JacobianPre"; 3817 DM dmAux = NULL; 3818 PetscDS prob, probAux = NULL; 3819 PetscSection sectionAux = NULL; 3820 Vec A; 3821 DMField coordField; 3822 PetscFEGeom *cgeomFEM; 3823 PetscQuadrature qGeom = NULL; 3824 Mat J = Jac, JP = JacP; 3825 PetscScalar *work, *u = NULL, *u_t = NULL, *a = NULL, *elemMat = NULL, *elemMatP = NULL, *elemMatD = NULL; 3826 PetscBool hasJac, hasPrec, hasDyn, assembleJac, isMatIS, isMatISP, *isFE, hasFV = PETSC_FALSE; 3827 const PetscInt *cells; 3828 PetscInt Nf, fieldI, fieldJ, maxDegree, numCells, cStart, cEnd, numChunks, chunkSize, chunk, totDim, totDimAux = 0, sz, wsz, off = 0, offCell = 0; 3829 PetscErrorCode ierr; 3830 3831 PetscFunctionBegin; 3832 CHKMEMQ; 3833 ierr = ISGetLocalSize(cellIS, &numCells);CHKERRQ(ierr); 3834 ierr = ISGetPointRange(cellIS, &cStart, &cEnd, &cells);CHKERRQ(ierr); 3835 ierr = PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr); 3836 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 3837 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 3838 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 3839 if (dmAux) { 3840 ierr = DMGetLocalSection(dmAux, §ionAux);CHKERRQ(ierr); 3841 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 3842 } 3843 /* Get flags */ 3844 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 3845 ierr = DMGetWorkArray(dm, Nf, MPIU_BOOL, &isFE);CHKERRQ(ierr); 3846 for (fieldI = 0; fieldI < Nf; ++fieldI) { 3847 PetscObject disc; 3848 PetscClassId id; 3849 ierr = PetscDSGetDiscretization(prob, fieldI, &disc);CHKERRQ(ierr); 3850 ierr = PetscObjectGetClassId(disc, &id);CHKERRQ(ierr); 3851 if (id == PETSCFE_CLASSID) {isFE[fieldI] = PETSC_TRUE;} 3852 else if (id == PETSCFV_CLASSID) {hasFV = PETSC_TRUE; isFE[fieldI] = PETSC_FALSE;} 3853 } 3854 ierr = PetscDSHasJacobian(prob, &hasJac);CHKERRQ(ierr); 3855 ierr = PetscDSHasJacobianPreconditioner(prob, &hasPrec);CHKERRQ(ierr); 3856 ierr = PetscDSHasDynamicJacobian(prob, &hasDyn);CHKERRQ(ierr); 3857 assembleJac = hasJac && hasPrec && (Jac != JacP) ? PETSC_TRUE : PETSC_FALSE; 3858 hasDyn = hasDyn && (X_tShift != 0.0) ? PETSC_TRUE : PETSC_FALSE; 3859 ierr = PetscObjectTypeCompare((PetscObject) Jac, MATIS, &isMatIS);CHKERRQ(ierr); 3860 ierr = PetscObjectTypeCompare((PetscObject) JacP, MATIS, &isMatISP);CHKERRQ(ierr); 3861 /* Setup input data and temp arrays (should be DMGetWorkArray) */ 3862 if (isMatISP || isMatISP) {ierr = DMPlexGetSubdomainSection(dm, &globalSection);CHKERRQ(ierr);} 3863 if (isMatIS) {ierr = MatISGetLocalMat(Jac, &J);CHKERRQ(ierr);} 3864 if (isMatISP) {ierr = MatISGetLocalMat(JacP, &JP);CHKERRQ(ierr);} 3865 if (hasFV) {ierr = MatSetOption(JP, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);CHKERRQ(ierr);} /* No allocated space for FV stuff, so ignore the zero entries */ 3866 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 3867 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 3868 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 3869 if (probAux) {ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr);} 3870 CHKMEMQ; 3871 /* Compute batch sizes */ 3872 if (isFE[0]) { 3873 PetscFE fe; 3874 PetscQuadrature q; 3875 PetscInt numQuadPoints, numBatches, batchSize, numBlocks, blockSize, Nb; 3876 3877 ierr = PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);CHKERRQ(ierr); 3878 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 3879 ierr = PetscQuadratureGetData(q, NULL, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 3880 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 3881 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 3882 blockSize = Nb*numQuadPoints; 3883 batchSize = numBlocks * blockSize; 3884 chunkSize = numBatches * batchSize; 3885 numChunks = numCells / chunkSize + numCells % chunkSize; 3886 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 3887 } else { 3888 chunkSize = numCells; 3889 numChunks = 1; 3890 } 3891 /* Get work space */ 3892 wsz = (((X?1:0) + (X_t?1:0) + (dmAux?1:0))*totDim + ((hasJac?1:0) + (hasPrec?1:0) + (hasDyn?1:0))*totDim*totDim)*chunkSize; 3893 ierr = DMGetWorkArray(dm, wsz, MPIU_SCALAR, &work);CHKERRQ(ierr); 3894 ierr = PetscArrayzero(work, wsz);CHKERRQ(ierr); 3895 off = 0; 3896 u = X ? (sz = chunkSize*totDim, off += sz, work+off-sz) : NULL; 3897 u_t = X_t ? (sz = chunkSize*totDim, off += sz, work+off-sz) : NULL; 3898 a = dmAux ? (sz = chunkSize*totDimAux, off += sz, work+off-sz) : NULL; 3899 elemMat = hasJac ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL; 3900 elemMatP = hasPrec ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL; 3901 elemMatD = hasDyn ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL; 3902 if (off != wsz) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error is workspace size %D should be %D", off, wsz); 3903 /* Setup geometry */ 3904 ierr = DMGetCoordinateField(dm, &coordField);CHKERRQ(ierr); 3905 ierr = DMFieldGetDegree(coordField, cellIS, NULL, &maxDegree);CHKERRQ(ierr); 3906 if (maxDegree <= 1) {ierr = DMFieldCreateDefaultQuadrature(coordField, cellIS, &qGeom);CHKERRQ(ierr);} 3907 if (!qGeom) { 3908 PetscFE fe; 3909 3910 ierr = PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);CHKERRQ(ierr); 3911 ierr = PetscFEGetQuadrature(fe, &qGeom);CHKERRQ(ierr); 3912 ierr = PetscObjectReference((PetscObject) qGeom);CHKERRQ(ierr); 3913 } 3914 ierr = DMSNESGetFEGeom(coordField, cellIS, qGeom, PETSC_FALSE, &cgeomFEM);CHKERRQ(ierr); 3915 /* Compute volume integrals */ 3916 if (assembleJac) {ierr = MatZeroEntries(J);CHKERRQ(ierr);} 3917 ierr = MatZeroEntries(JP);CHKERRQ(ierr); 3918 for (chunk = 0; chunk < numChunks; ++chunk, offCell += chunkSize) { 3919 const PetscInt Ncell = PetscMin(chunkSize, numCells - offCell); 3920 PetscInt c; 3921 3922 /* Extract values */ 3923 for (c = 0; c < Ncell; ++c) { 3924 const PetscInt cell = cells ? cells[c+offCell] : c+offCell; 3925 PetscScalar *x = NULL, *x_t = NULL; 3926 PetscInt i; 3927 3928 if (X) { 3929 ierr = DMPlexVecGetClosure(dm, section, X, cell, NULL, &x);CHKERRQ(ierr); 3930 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 3931 ierr = DMPlexVecRestoreClosure(dm, section, X, cell, NULL, &x);CHKERRQ(ierr); 3932 } 3933 if (X_t) { 3934 ierr = DMPlexVecGetClosure(dm, section, X_t, cell, NULL, &x_t);CHKERRQ(ierr); 3935 for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i]; 3936 ierr = DMPlexVecRestoreClosure(dm, section, X_t, cell, NULL, &x_t);CHKERRQ(ierr); 3937 } 3938 if (dmAux) { 3939 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, cell, NULL, &x);CHKERRQ(ierr); 3940 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 3941 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, cell, NULL, &x);CHKERRQ(ierr); 3942 } 3943 } 3944 CHKMEMQ; 3945 for (fieldI = 0; fieldI < Nf; ++fieldI) { 3946 PetscFE fe; 3947 ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr); 3948 for (fieldJ = 0; fieldJ < Nf; ++fieldJ) { 3949 if (hasJac) {ierr = PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMat);CHKERRQ(ierr);} 3950 if (hasPrec) {ierr = PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_PRE, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMatP);CHKERRQ(ierr);} 3951 if (hasDyn) {ierr = PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMatD);CHKERRQ(ierr);} 3952 } 3953 /* For finite volume, add the identity */ 3954 if (!isFE[fieldI]) { 3955 PetscFV fv; 3956 PetscInt eOffset = 0, Nc, fc, foff; 3957 3958 ierr = PetscDSGetFieldOffset(prob, fieldI, &foff);CHKERRQ(ierr); 3959 ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fv);CHKERRQ(ierr); 3960 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 3961 for (c = 0; c < chunkSize; ++c, eOffset += totDim*totDim) { 3962 for (fc = 0; fc < Nc; ++fc) { 3963 const PetscInt i = foff + fc; 3964 if (hasJac) {elemMat [eOffset+i*totDim+i] = 1.0;} 3965 if (hasPrec) {elemMatP[eOffset+i*totDim+i] = 1.0;} 3966 } 3967 } 3968 } 3969 } 3970 CHKMEMQ; 3971 /* Add contribution from X_t */ 3972 if (hasDyn) {for (c = 0; c < chunkSize*totDim*totDim; ++c) elemMat[c] += X_tShift*elemMatD[c];} 3973 /* Insert values into matrix */ 3974 for (c = 0; c < Ncell; ++c) { 3975 const PetscInt cell = cells ? cells[c+offCell] : c+offCell; 3976 if (mesh->printFEM > 1) { 3977 if (hasJac) {ierr = DMPrintCellMatrix(cell, name, totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);CHKERRQ(ierr);} 3978 if (hasPrec) {ierr = DMPrintCellMatrix(cell, nameP, totDim, totDim, &elemMatP[(c-cStart)*totDim*totDim]);CHKERRQ(ierr);} 3979 } 3980 if (assembleJac) {ierr = DMPlexMatSetClosure(dm, section, globalSection, Jac, cell, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);CHKERRQ(ierr);} 3981 ierr = DMPlexMatSetClosure(dm, section, globalSection, JP, cell, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);CHKERRQ(ierr); 3982 } 3983 CHKMEMQ; 3984 } 3985 /* Cleanup */ 3986 ierr = DMSNESRestoreFEGeom(coordField, cellIS, qGeom, PETSC_FALSE, &cgeomFEM);CHKERRQ(ierr); 3987 ierr = PetscQuadratureDestroy(&qGeom);CHKERRQ(ierr); 3988 if (hasFV) {ierr = MatSetOption(JacP, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);CHKERRQ(ierr);} 3989 ierr = DMRestoreWorkArray(dm, Nf, MPIU_BOOL, &isFE);CHKERRQ(ierr); 3990 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); 3991 /* Compute boundary integrals */ 3992 /* ierr = DMPlexComputeBdJacobian_Internal(dm, X, X_t, t, X_tShift, Jac, JacP, ctx);CHKERRQ(ierr); */ 3993 /* Assemble matrix */ 3994 if (assembleJac) {ierr = MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);ierr = MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);} 3995 ierr = MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);ierr = MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3996 ierr = PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr); 3997 CHKMEMQ; 3998 PetscFunctionReturn(0); 3999 } 4000