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