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