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