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