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