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