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