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