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