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