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