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