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