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