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