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