Actual source code: plexfem.c
petsc-3.14.6 2021-03-30
1: #include <petsc/private/dmpleximpl.h>
2: #include <petscsf.h>
4: #include <petsc/private/hashsetij.h>
5: #include <petsc/private/petscfeimpl.h>
6: #include <petsc/private/petscfvimpl.h>
8: static PetscErrorCode DMPlexConvertPlex(DM dm, DM *plex, PetscBool copy)
9: {
10: PetscBool isPlex;
14: PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);
15: if (isPlex) {
16: *plex = dm;
17: PetscObjectReference((PetscObject) dm);
18: } else {
19: PetscObjectQuery((PetscObject) dm, "dm_plex", (PetscObject *) plex);
20: if (!*plex) {
21: DMConvert(dm,DMPLEX,plex);
22: PetscObjectCompose((PetscObject) dm, "dm_plex", (PetscObject) *plex);
23: if (copy) {
24: const char *comps[3] = {"A", "dmAux"};
25: PetscObject obj;
26: PetscInt i;
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) {(*link->ddhook)(dm, *plex, link->ctx);}
33: }
34: }
35: for (i = 0; i < 3; i++) {
36: PetscObjectQuery((PetscObject) dm, comps[i], &obj);
37: PetscObjectCompose((PetscObject) *plex, comps[i], obj);
38: }
39: }
40: } else {
41: PetscObjectReference((PetscObject) *plex);
42: }
43: }
44: return(0);
45: }
47: static PetscErrorCode PetscContainerUserDestroy_PetscFEGeom (void *ctx)
48: {
49: PetscFEGeom *geom = (PetscFEGeom *) ctx;
53: PetscFEGeomDestroy(&geom);
54: return(0);
55: }
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;
65: PetscObjectGetId((PetscObject)quad,&id);
66: PetscSNPrintf(composeStr, 32, "DMPlexGetFEGeom_%x\n", id);
67: PetscObjectQuery((PetscObject) pointIS, composeStr, (PetscObject *) &container);
68: if (container) {
69: PetscContainerGetPointer(container, (void **) geom);
70: } else {
71: DMFieldCreateFEGeom(coordField, pointIS, quad, faceData, geom);
72: PetscContainerCreate(PETSC_COMM_SELF,&container);
73: PetscContainerSetPointer(container, (void *) *geom);
74: PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom);
75: PetscObjectCompose((PetscObject) pointIS, composeStr, (PetscObject) container);
76: PetscContainerDestroy(&container);
77: }
78: return(0);
79: }
81: static PetscErrorCode DMPlexRestoreFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
82: {
84: *geom = NULL;
85: return(0);
86: }
88: /*@
89: DMPlexGetScale - Get the scale for the specified fundamental unit
91: Not collective
93: Input Arguments:
94: + dm - the DM
95: - unit - The SI unit
97: Output Argument:
98: . scale - The value used to scale all quantities with this unit
100: Level: advanced
102: .seealso: DMPlexSetScale(), PetscUnit
103: @*/
104: PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale)
105: {
106: DM_Plex *mesh = (DM_Plex*) dm->data;
111: *scale = mesh->scale[unit];
112: return(0);
113: }
115: /*@
116: DMPlexSetScale - Set the scale for the specified fundamental unit
118: Not collective
120: Input Arguments:
121: + dm - the DM
122: . unit - The SI unit
123: - scale - The value used to scale all quantities with this unit
125: Level: advanced
127: .seealso: DMPlexGetScale(), PetscUnit
128: @*/
129: PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale)
130: {
131: DM_Plex *mesh = (DM_Plex*) dm->data;
135: mesh->scale[unit] = scale;
136: return(0);
137: }
139: static PetscErrorCode DMPlexProjectRigidBody_Private(PetscInt dim, PetscReal t, const PetscReal X[], PetscInt Nc, PetscScalar *mode, void *ctx)
140: {
141: const PetscInt eps[3][3][3] = {{{0, 0, 0}, {0, 0, 1}, {0, -1, 0}}, {{0, 0, -1}, {0, 0, 0}, {1, 0, 0}}, {{0, 1, 0}, {-1, 0, 0}, {0, 0, 0}}};
142: PetscInt *ctxInt = (PetscInt *) ctx;
143: PetscInt dim2 = ctxInt[0];
144: PetscInt d = ctxInt[1];
145: PetscInt i, j, k = dim > 2 ? d - dim : d;
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: return(0);
160: }
162: /*@
163: DMPlexCreateRigidBody - For the default global section, create rigid body modes by function space interpolation
165: Collective on dm
167: Input Arguments:
168: + dm - the DM
169: - field - The field number for the rigid body space, or 0 for the default
171: Output Argument:
172: . sp - the null space
174: Note: This is necessary to provide a suitable coarse space for algebraic multigrid
176: Level: advanced
178: .seealso: MatNullSpaceCreate(), PCGAMG
179: @*/
180: PetscErrorCode DMPlexCreateRigidBody(DM dm, PetscInt field, MatNullSpace *sp)
181: {
182: PetscErrorCode (**func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *);
183: MPI_Comm comm;
184: Vec mode[6];
185: PetscSection section, globalSection;
186: PetscInt dim, dimEmbed, Nf, n, m, mmin, d, i, j;
187: PetscErrorCode ierr;
190: PetscObjectGetComm((PetscObject) dm, &comm);
191: DMGetDimension(dm, &dim);
192: DMGetCoordinateDim(dm, &dimEmbed);
193: DMGetNumFields(dm, &Nf);
194: if (Nf && (field < 0 || field >= Nf)) SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Field %D is not in [0, Nf)", field, Nf);
195: if (dim == 1 && Nf < 2) {
196: MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);
197: return(0);
198: }
199: DMGetLocalSection(dm, §ion);
200: DMGetGlobalSection(dm, &globalSection);
201: PetscSectionGetConstrainedStorageSize(globalSection, &n);
202: PetscCalloc1(Nf, &func);
203: m = (dim*(dim+1))/2;
204: VecCreate(comm, &mode[0]);
205: VecSetSizes(mode[0], n, PETSC_DETERMINE);
206: VecSetUp(mode[0]);
207: VecGetSize(mode[0], &n);
208: mmin = PetscMin(m, n);
209: func[field] = DMPlexProjectRigidBody_Private;
210: for (i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
211: for (d = 0; d < m; d++) {
212: PetscInt ctx[2];
213: void *voidctx = (void *) (&ctx[0]);
215: ctx[0] = dimEmbed;
216: ctx[1] = d;
217: DMProjectFunction(dm, 0.0, func, &voidctx, INSERT_VALUES, mode[d]);
218: }
219: /* Orthonormalize system */
220: for (i = 0; i < mmin; ++i) {
221: PetscScalar dots[6];
223: VecNormalize(mode[i], NULL);
224: VecMDot(mode[i], mmin-i-1, mode+i+1, dots+i+1);
225: for (j = i+1; j < mmin; ++j) {
226: dots[j] *= -1.0;
227: VecAXPY(mode[j], dots[j], mode[i]);
228: }
229: }
230: MatNullSpaceCreate(comm, PETSC_FALSE, mmin, mode, sp);
231: for (i = 0; i < m; ++i) {VecDestroy(&mode[i]);}
232: PetscFree(func);
233: return(0);
234: }
236: /*@
237: DMPlexCreateRigidBodies - For the default global section, create rigid body modes by function space interpolation
239: Collective on dm
241: Input Arguments:
242: + dm - the DM
243: . nb - The number of bodies
244: . label - The DMLabel marking each domain
245: . nids - The number of ids per body
246: - ids - An array of the label ids in sequence for each domain
248: Output Argument:
249: . sp - the null space
251: Note: This is necessary to provide a suitable coarse space for algebraic multigrid
253: Level: advanced
255: .seealso: MatNullSpaceCreate()
256: @*/
257: PetscErrorCode DMPlexCreateRigidBodies(DM dm, PetscInt nb, DMLabel label, const PetscInt nids[], const PetscInt ids[], MatNullSpace *sp)
258: {
259: MPI_Comm comm;
260: PetscSection section, globalSection;
261: Vec *mode;
262: PetscScalar *dots;
263: PetscInt dim, dimEmbed, n, m, b, d, i, j, off;
267: PetscObjectGetComm((PetscObject)dm,&comm);
268: DMGetDimension(dm, &dim);
269: DMGetCoordinateDim(dm, &dimEmbed);
270: DMGetLocalSection(dm, §ion);
271: DMGetGlobalSection(dm, &globalSection);
272: PetscSectionGetConstrainedStorageSize(globalSection, &n);
273: m = nb * (dim*(dim+1))/2;
274: PetscMalloc2(m, &mode, m, &dots);
275: VecCreate(comm, &mode[0]);
276: VecSetSizes(mode[0], n, PETSC_DETERMINE);
277: VecSetUp(mode[0]);
278: for (i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
279: for (b = 0, off = 0; b < nb; ++b) {
280: for (d = 0; d < m/nb; ++d) {
281: PetscInt ctx[2];
282: PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private;
283: void *voidctx = (void *) (&ctx[0]);
285: ctx[0] = dimEmbed;
286: ctx[1] = d;
287: DMProjectFunctionLabel(dm, 0.0, label, nids[b], &ids[off], 0, NULL, &func, &voidctx, INSERT_VALUES, mode[d]);
288: off += nids[b];
289: }
290: }
291: /* Orthonormalize system */
292: for (i = 0; i < m; ++i) {
293: PetscScalar dots[6];
295: VecNormalize(mode[i], NULL);
296: VecMDot(mode[i], m-i-1, mode+i+1, dots+i+1);
297: for (j = i+1; j < m; ++j) {
298: dots[j] *= -1.0;
299: VecAXPY(mode[j], dots[j], mode[i]);
300: }
301: }
302: MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);
303: for (i = 0; i< m; ++i) {VecDestroy(&mode[i]);}
304: PetscFree2(mode, dots);
305: return(0);
306: }
308: /*@
309: DMPlexSetMaxProjectionHeight - In DMPlexProjectXXXLocal() functions, the projected values of a basis function's dofs
310: are computed by associating the basis function with one of the mesh points in its transitively-closed support, and
311: evaluating the dual space basis of that point. A basis function is associated with the point in its
312: transitively-closed support whose mesh height is highest (w.r.t. DAG height), but not greater than the maximum
313: projection height, which is set with this function. By default, the maximum projection height is zero, which means
314: that only mesh cells are used to project basis functions. A height of one, for example, evaluates a cell-interior
315: basis functions using its cells dual space basis, but all other basis functions with the dual space basis of a face.
317: Input Parameters:
318: + dm - the DMPlex object
319: - height - the maximum projection height >= 0
321: Level: advanced
323: .seealso: DMPlexGetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
324: @*/
325: PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height)
326: {
327: DM_Plex *plex = (DM_Plex *) dm->data;
331: plex->maxProjectionHeight = height;
332: return(0);
333: }
335: /*@
336: DMPlexGetMaxProjectionHeight - Get the maximum height (w.r.t. DAG) of mesh points used to evaluate dual bases in
337: DMPlexProjectXXXLocal() functions.
339: Input Parameters:
340: . dm - the DMPlex object
342: Output Parameters:
343: . height - the maximum projection height
345: Level: intermediate
347: .seealso: DMPlexSetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
348: @*/
349: PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height)
350: {
351: DM_Plex *plex = (DM_Plex *) dm->data;
355: *height = plex->maxProjectionHeight;
356: return(0);
357: }
359: typedef struct {
360: PetscReal alpha; /* The first Euler angle, and in 2D the only one */
361: PetscReal beta; /* The second Euler angle */
362: PetscReal gamma; /* The third Euler angle */
363: PetscInt dim; /* The dimension of R */
364: PetscScalar *R; /* The rotation matrix, transforming a vector in the local basis to the global basis */
365: PetscScalar *RT; /* The transposed rotation matrix, transforming a vector in the global basis to the local basis */
366: } RotCtx;
368: /*
369: Note: Following https://en.wikipedia.org/wiki/Euler_angles, we will specify Euler angles by extrinsic rotations, meaning that
370: 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:
371: $ The XYZ system rotates about the z axis by alpha. The X axis is now at angle alpha with respect to the x axis.
372: $ 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.
373: $ The XYZ system rotates a third time about the z axis by gamma.
374: */
375: static PetscErrorCode DMPlexBasisTransformSetUp_Rotation_Internal(DM dm, void *ctx)
376: {
377: RotCtx *rc = (RotCtx *) ctx;
378: PetscInt dim = rc->dim;
379: PetscReal c1, s1, c2, s2, c3, s3;
383: PetscMalloc2(PetscSqr(dim), &rc->R, PetscSqr(dim), &rc->RT);
384: switch (dim) {
385: case 2:
386: c1 = PetscCosReal(rc->alpha);s1 = PetscSinReal(rc->alpha);
387: rc->R[0] = c1;rc->R[1] = s1;
388: rc->R[2] = -s1;rc->R[3] = c1;
389: PetscArraycpy(rc->RT, rc->R, PetscSqr(dim));
390: DMPlex_Transpose2D_Internal(rc->RT);
391: break;
392: case 3:
393: c1 = PetscCosReal(rc->alpha);s1 = PetscSinReal(rc->alpha);
394: c2 = PetscCosReal(rc->beta); s2 = PetscSinReal(rc->beta);
395: c3 = PetscCosReal(rc->gamma);s3 = PetscSinReal(rc->gamma);
396: rc->R[0] = c1*c3 - c2*s1*s3;rc->R[1] = c3*s1 + c1*c2*s3;rc->R[2] = s2*s3;
397: rc->R[3] = -c1*s3 - c2*c3*s1;rc->R[4] = c1*c2*c3 - s1*s3; rc->R[5] = c3*s2;
398: rc->R[6] = s1*s2; rc->R[7] = -c1*s2; rc->R[8] = c2;
399: PetscArraycpy(rc->RT, rc->R, PetscSqr(dim));
400: DMPlex_Transpose3D_Internal(rc->RT);
401: break;
402: default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not supported", dim);
403: }
404: return(0);
405: }
407: static PetscErrorCode DMPlexBasisTransformDestroy_Rotation_Internal(DM dm, void *ctx)
408: {
409: RotCtx *rc = (RotCtx *) ctx;
413: PetscFree2(rc->R, rc->RT);
414: PetscFree(rc);
415: return(0);
416: }
418: static PetscErrorCode DMPlexBasisTransformGetMatrix_Rotation_Internal(DM dm, const PetscReal x[], PetscBool l2g, const PetscScalar **A, void *ctx)
419: {
420: RotCtx *rc = (RotCtx *) ctx;
424: if (l2g) {*A = rc->R;}
425: else {*A = rc->RT;}
426: return(0);
427: }
429: PetscErrorCode DMPlexBasisTransformApplyReal_Internal(DM dm, const PetscReal x[], PetscBool l2g, PetscInt dim, const PetscReal *y, PetscReal *z, void *ctx)
430: {
434: #if defined(PETSC_USE_COMPLEX)
435: switch (dim) {
436: case 2:
437: {
438: PetscScalar yt[2], zt[2];
440: yt[0] = y[0]; yt[1] = y[1];
441: DMPlexBasisTransformApply_Internal(dm, x, l2g, dim, yt, zt, ctx);
442: z[0] = PetscRealPart(zt[0]); z[1] = PetscRealPart(zt[1]);
443: }
444: break;
445: case 3:
446: {
447: PetscScalar yt[3], zt[3];
449: yt[0] = y[0]; yt[1] = y[1]; yt[2] = y[2];
450: DMPlexBasisTransformApply_Internal(dm, x, l2g, dim, yt, zt, ctx);
451: z[0] = PetscRealPart(zt[0]); z[1] = PetscRealPart(zt[1]); z[2] = PetscRealPart(zt[2]);
452: }
453: break;
454: }
455: #else
456: DMPlexBasisTransformApply_Internal(dm, x, l2g, dim, y, z, ctx);
457: #endif
458: return(0);
459: }
461: PetscErrorCode DMPlexBasisTransformApply_Internal(DM dm, const PetscReal x[], PetscBool l2g, PetscInt dim, const PetscScalar *y, PetscScalar *z, void *ctx)
462: {
463: const PetscScalar *A;
464: PetscErrorCode ierr;
467: (*dm->transformGetMatrix)(dm, x, l2g, &A, ctx);
468: switch (dim) {
469: case 2: DMPlex_Mult2D_Internal(A, 1, y, z);break;
470: case 3: DMPlex_Mult3D_Internal(A, 1, y, z);break;
471: }
472: return(0);
473: }
475: static PetscErrorCode DMPlexBasisTransformField_Internal(DM dm, DM tdm, Vec tv, PetscInt p, PetscInt f, PetscBool l2g, PetscScalar *a)
476: {
477: PetscSection ts;
478: const PetscScalar *ta, *tva;
479: PetscInt dof;
480: PetscErrorCode ierr;
483: DMGetLocalSection(tdm, &ts);
484: PetscSectionGetFieldDof(ts, p, f, &dof);
485: VecGetArrayRead(tv, &ta);
486: DMPlexPointLocalFieldRead(tdm, p, f, ta, (void *) &tva);
487: if (l2g) {
488: switch (dof) {
489: case 4: DMPlex_Mult2D_Internal(tva, 1, a, a);break;
490: case 9: DMPlex_Mult3D_Internal(tva, 1, a, a);break;
491: }
492: } else {
493: switch (dof) {
494: case 4: DMPlex_MultTranspose2D_Internal(tva, 1, a, a);break;
495: case 9: DMPlex_MultTranspose3D_Internal(tva, 1, a, a);break;
496: }
497: }
498: VecRestoreArrayRead(tv, &ta);
499: return(0);
500: }
502: static PetscErrorCode DMPlexBasisTransformFieldTensor_Internal(DM dm, DM tdm, Vec tv, PetscInt pf, PetscInt f, PetscInt pg, PetscInt g, PetscBool l2g, PetscInt lda, PetscScalar *a)
503: {
504: PetscSection s, ts;
505: const PetscScalar *ta, *tvaf, *tvag;
506: PetscInt fdof, gdof, fpdof, gpdof;
507: PetscErrorCode ierr;
510: DMGetLocalSection(dm, &s);
511: DMGetLocalSection(tdm, &ts);
512: PetscSectionGetFieldDof(s, pf, f, &fpdof);
513: PetscSectionGetFieldDof(s, pg, g, &gpdof);
514: PetscSectionGetFieldDof(ts, pf, f, &fdof);
515: PetscSectionGetFieldDof(ts, pg, g, &gdof);
516: VecGetArrayRead(tv, &ta);
517: DMPlexPointLocalFieldRead(tdm, pf, f, ta, (void *) &tvaf);
518: DMPlexPointLocalFieldRead(tdm, pg, g, ta, (void *) &tvag);
519: if (l2g) {
520: switch (fdof) {
521: case 4: DMPlex_MatMult2D_Internal(tvaf, gpdof, lda, a, a);break;
522: case 9: DMPlex_MatMult3D_Internal(tvaf, gpdof, lda, a, a);break;
523: }
524: switch (gdof) {
525: case 4: DMPlex_MatMultTransposeLeft2D_Internal(tvag, fpdof, lda, a, a);break;
526: case 9: DMPlex_MatMultTransposeLeft3D_Internal(tvag, fpdof, lda, a, a);break;
527: }
528: } else {
529: switch (fdof) {
530: case 4: DMPlex_MatMultTranspose2D_Internal(tvaf, gpdof, lda, a, a);break;
531: case 9: DMPlex_MatMultTranspose3D_Internal(tvaf, gpdof, lda, a, a);break;
532: }
533: switch (gdof) {
534: case 4: DMPlex_MatMultLeft2D_Internal(tvag, fpdof, lda, a, a);break;
535: case 9: DMPlex_MatMultLeft3D_Internal(tvag, fpdof, lda, a, a);break;
536: }
537: }
538: VecRestoreArrayRead(tv, &ta);
539: return(0);
540: }
542: PetscErrorCode DMPlexBasisTransformPoint_Internal(DM dm, DM tdm, Vec tv, PetscInt p, PetscBool fieldActive[], PetscBool l2g, PetscScalar *a)
543: {
544: PetscSection s;
545: PetscSection clSection;
546: IS clPoints;
547: const PetscInt *clp;
548: PetscInt *points = NULL;
549: PetscInt Nf, f, Np, cp, dof, d = 0;
550: PetscErrorCode ierr;
553: DMGetLocalSection(dm, &s);
554: PetscSectionGetNumFields(s, &Nf);
555: DMPlexGetCompressedClosure(dm, s, p, &Np, &points, &clSection, &clPoints, &clp);
556: for (f = 0; f < Nf; ++f) {
557: for (cp = 0; cp < Np*2; cp += 2) {
558: PetscSectionGetFieldDof(s, points[cp], f, &dof);
559: if (!dof) continue;
560: if (fieldActive[f]) {DMPlexBasisTransformField_Internal(dm, tdm, tv, points[cp], f, l2g, &a[d]);}
561: d += dof;
562: }
563: }
564: DMPlexRestoreCompressedClosure(dm, s, p, &Np, &points, &clSection, &clPoints, &clp);
565: return(0);
566: }
568: PetscErrorCode DMPlexBasisTransformPointTensor_Internal(DM dm, DM tdm, Vec tv, PetscInt p, PetscBool l2g, PetscInt lda, PetscScalar *a)
569: {
570: PetscSection s;
571: PetscSection clSection;
572: IS clPoints;
573: const PetscInt *clp;
574: PetscInt *points = NULL;
575: PetscInt Nf, f, g, Np, cpf, cpg, fdof, gdof, r, c = 0;
576: PetscErrorCode ierr;
579: DMGetLocalSection(dm, &s);
580: PetscSectionGetNumFields(s, &Nf);
581: DMPlexGetCompressedClosure(dm, s, p, &Np, &points, &clSection, &clPoints, &clp);
582: for (f = 0, r = 0; f < Nf; ++f) {
583: for (cpf = 0; cpf < Np*2; cpf += 2) {
584: PetscSectionGetFieldDof(s, points[cpf], f, &fdof);
585: for (g = 0, c = 0; g < Nf; ++g) {
586: for (cpg = 0; cpg < Np*2; cpg += 2) {
587: PetscSectionGetFieldDof(s, points[cpg], g, &gdof);
588: DMPlexBasisTransformFieldTensor_Internal(dm, tdm, tv, points[cpf], f, points[cpg], g, l2g, lda, &a[r*lda+c]);
589: c += gdof;
590: }
591: }
592: if (c != lda) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of columns %D should be %D", c, lda);
593: r += fdof;
594: }
595: }
596: if (r != lda) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of rows %D should be %D", c, lda);
597: DMPlexRestoreCompressedClosure(dm, s, p, &Np, &points, &clSection, &clPoints, &clp);
598: return(0);
599: }
601: static PetscErrorCode DMPlexBasisTransform_Internal(DM dm, Vec lv, PetscBool l2g)
602: {
603: DM tdm;
604: Vec tv;
605: PetscSection ts, s;
606: const PetscScalar *ta;
607: PetscScalar *a, *va;
608: PetscInt pStart, pEnd, p, Nf, f;
609: PetscErrorCode ierr;
612: DMGetBasisTransformDM_Internal(dm, &tdm);
613: DMGetBasisTransformVec_Internal(dm, &tv);
614: DMGetLocalSection(tdm, &ts);
615: DMGetLocalSection(dm, &s);
616: PetscSectionGetChart(s, &pStart, &pEnd);
617: PetscSectionGetNumFields(s, &Nf);
618: VecGetArray(lv, &a);
619: VecGetArrayRead(tv, &ta);
620: for (p = pStart; p < pEnd; ++p) {
621: for (f = 0; f < Nf; ++f) {
622: DMPlexPointLocalFieldRef(dm, p, f, a, (void *) &va);
623: DMPlexBasisTransformField_Internal(dm, tdm, tv, p, f, l2g, va);
624: }
625: }
626: VecRestoreArray(lv, &a);
627: VecRestoreArrayRead(tv, &ta);
628: return(0);
629: }
631: /*@
632: DMPlexGlobalToLocalBasis - Transform the values in the given local vector from the global basis to the local basis
634: Input Parameters:
635: + dm - The DM
636: - lv - A local vector with values in the global basis
638: Output Parameters:
639: . lv - A local vector with values in the local basis
641: 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.
643: Level: developer
645: .seealso: DMPlexLocalToGlobalBasis(), DMGetLocalSection(), DMPlexCreateBasisRotation()
646: @*/
647: PetscErrorCode DMPlexGlobalToLocalBasis(DM dm, Vec lv)
648: {
654: DMPlexBasisTransform_Internal(dm, lv, PETSC_FALSE);
655: return(0);
656: }
658: /*@
659: DMPlexLocalToGlobalBasis - Transform the values in the given local vector from the local basis to the global basis
661: Input Parameters:
662: + dm - The DM
663: - lv - A local vector with values in the local basis
665: Output Parameters:
666: . lv - A local vector with values in the global basis
668: 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.
670: Level: developer
672: .seealso: DMPlexGlobalToLocalBasis(), DMGetLocalSection(), DMPlexCreateBasisRotation()
673: @*/
674: PetscErrorCode DMPlexLocalToGlobalBasis(DM dm, Vec lv)
675: {
681: DMPlexBasisTransform_Internal(dm, lv, PETSC_TRUE);
682: return(0);
683: }
685: /*@
686: DMPlexCreateBasisRotation - Create an internal transformation from the global basis, used to specify boundary conditions
687: and global solutions, to a local basis, appropriate for discretization integrals and assembly.
689: Input Parameters:
690: + dm - The DM
691: . alpha - The first Euler angle, and in 2D the only one
692: . beta - The second Euler angle
693: - gamma - The third Euler angle
695: Note: Following https://en.wikipedia.org/wiki/Euler_angles, we will specify Euler angles by extrinsic rotations, meaning that
696: 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:
697: $ The XYZ system rotates about the z axis by alpha. The X axis is now at angle alpha with respect to the x axis.
698: $ 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.
699: $ The XYZ system rotates a third time about the z axis by gamma.
701: Level: developer
703: .seealso: DMPlexGlobalToLocalBasis(), DMPlexLocalToGlobalBasis()
704: @*/
705: PetscErrorCode DMPlexCreateBasisRotation(DM dm, PetscReal alpha, PetscReal beta, PetscReal gamma)
706: {
707: RotCtx *rc;
708: PetscInt cdim;
711: DMGetCoordinateDim(dm, &cdim);
712: PetscMalloc1(1, &rc);
713: dm->transformCtx = rc;
714: dm->transformSetUp = DMPlexBasisTransformSetUp_Rotation_Internal;
715: dm->transformDestroy = DMPlexBasisTransformDestroy_Rotation_Internal;
716: dm->transformGetMatrix = DMPlexBasisTransformGetMatrix_Rotation_Internal;
717: rc->dim = cdim;
718: rc->alpha = alpha;
719: rc->beta = beta;
720: rc->gamma = gamma;
721: (*dm->transformSetUp)(dm, dm->transformCtx);
722: DMConstructBasisTransform_Internal(dm);
723: return(0);
724: }
726: /*@C
727: DMPlexInsertBoundaryValuesEssential - Insert boundary values into a local vector using a function of the coordinates
729: Input Parameters:
730: + dm - The DM, with a PetscDS that matches the problem being constrained
731: . time - The time
732: . field - The field to constrain
733: . Nc - The number of constrained field components, or 0 for all components
734: . comps - An array of constrained component numbers, or NULL for all components
735: . label - The DMLabel defining constrained points
736: . numids - The number of DMLabel ids for constrained points
737: . ids - An array of ids for constrained points
738: . func - A pointwise function giving boundary values
739: - ctx - An optional user context for bcFunc
741: Output Parameter:
742: . locX - A local vector to receives the boundary values
744: Level: developer
746: .seealso: DMPlexInsertBoundaryValuesEssentialField(), DMPlexInsertBoundaryValuesEssentialBdField(), DMAddBoundary()
747: @*/
748: 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)
749: {
750: PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
751: void **ctxs;
752: PetscInt numFields;
753: PetscErrorCode ierr;
756: DMGetNumFields(dm, &numFields);
757: PetscCalloc2(numFields,&funcs,numFields,&ctxs);
758: funcs[field] = func;
759: ctxs[field] = ctx;
760: DMProjectFunctionLabelLocal(dm, time, label, numids, ids, Nc, comps, funcs, ctxs, INSERT_BC_VALUES, locX);
761: PetscFree2(funcs,ctxs);
762: return(0);
763: }
765: /*@C
766: DMPlexInsertBoundaryValuesEssentialField - Insert boundary values into a local vector using a function of the coordinates and field data
768: Input Parameters:
769: + dm - The DM, with a PetscDS that matches the problem being constrained
770: . time - The time
771: . locU - A local vector with the input solution values
772: . field - The field to constrain
773: . Nc - The number of constrained field components, or 0 for all components
774: . comps - An array of constrained component numbers, or NULL for all components
775: . label - The DMLabel defining constrained points
776: . numids - The number of DMLabel ids for constrained points
777: . ids - An array of ids for constrained points
778: . func - A pointwise function giving boundary values
779: - ctx - An optional user context for bcFunc
781: Output Parameter:
782: . locX - A local vector to receives the boundary values
784: Level: developer
786: .seealso: DMPlexInsertBoundaryValuesEssential(), DMPlexInsertBoundaryValuesEssentialBdField(), DMAddBoundary()
787: @*/
788: PetscErrorCode DMPlexInsertBoundaryValuesEssentialField(DM dm, PetscReal time, Vec locU, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[],
789: void (*func)(PetscInt, PetscInt, PetscInt,
790: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
791: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
792: PetscReal, const PetscReal[], PetscInt, const PetscScalar[],
793: PetscScalar[]),
794: void *ctx, Vec locX)
795: {
796: void (**funcs)(PetscInt, PetscInt, PetscInt,
797: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
798: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
799: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
800: void **ctxs;
801: PetscInt numFields;
802: PetscErrorCode ierr;
805: DMGetNumFields(dm, &numFields);
806: PetscCalloc2(numFields,&funcs,numFields,&ctxs);
807: funcs[field] = func;
808: ctxs[field] = ctx;
809: DMProjectFieldLabelLocal(dm, time, label, numids, ids, Nc, comps, locU, funcs, INSERT_BC_VALUES, locX);
810: PetscFree2(funcs,ctxs);
811: return(0);
812: }
814: /*@C
815: DMPlexInsertBoundaryValuesEssentialBdField - Insert boundary values into a local vector using a function of the coodinates and boundary field data
817: Collective on dm
819: Input Parameters:
820: + dm - The DM, with a PetscDS that matches the problem being constrained
821: . time - The time
822: . locU - A local vector with the input solution values
823: . field - The field to constrain
824: . Nc - The number of constrained field components, or 0 for all components
825: . comps - An array of constrained component numbers, or NULL for all components
826: . label - The DMLabel defining constrained points
827: . numids - The number of DMLabel ids for constrained points
828: . ids - An array of ids for constrained points
829: . func - A pointwise function giving boundary values, the calling sequence is given in DMProjectBdFieldLabelLocal()
830: - ctx - An optional user context for bcFunc
832: Output Parameter:
833: . locX - A local vector to receive the boundary values
835: Level: developer
837: .seealso: DMProjectBdFieldLabelLocal(), DMPlexInsertBoundaryValuesEssential(), DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary()
838: @*/
839: PetscErrorCode DMPlexInsertBoundaryValuesEssentialBdField(DM dm, PetscReal time, Vec locU, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[],
840: void (*func)(PetscInt, PetscInt, PetscInt,
841: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
842: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
843: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[],
844: PetscScalar[]),
845: void *ctx, Vec locX)
846: {
847: void (**funcs)(PetscInt, PetscInt, PetscInt,
848: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
849: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
850: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
851: void **ctxs;
852: PetscInt numFields;
853: PetscErrorCode ierr;
856: DMGetNumFields(dm, &numFields);
857: PetscCalloc2(numFields,&funcs,numFields,&ctxs);
858: funcs[field] = func;
859: ctxs[field] = ctx;
860: DMProjectBdFieldLabelLocal(dm, time, label, numids, ids, Nc, comps, locU, funcs, INSERT_BC_VALUES, locX);
861: PetscFree2(funcs,ctxs);
862: return(0);
863: }
865: /*@C
866: DMPlexInsertBoundaryValuesRiemann - Insert boundary values into a local vector
868: Input Parameters:
869: + dm - The DM, with a PetscDS that matches the problem being constrained
870: . time - The time
871: . faceGeometry - A vector with the FVM face geometry information
872: . cellGeometry - A vector with the FVM cell geometry information
873: . Grad - A vector with the FVM cell gradient information
874: . field - The field to constrain
875: . Nc - The number of constrained field components, or 0 for all components
876: . comps - An array of constrained component numbers, or NULL for all components
877: . label - The DMLabel defining constrained points
878: . numids - The number of DMLabel ids for constrained points
879: . ids - An array of ids for constrained points
880: . func - A pointwise function giving boundary values
881: - ctx - An optional user context for bcFunc
883: Output Parameter:
884: . locX - A local vector to receives the boundary values
886: Note: This implementation currently ignores the numcomps/comps argument from DMAddBoundary()
888: Level: developer
890: .seealso: DMPlexInsertBoundaryValuesEssential(), DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary()
891: @*/
892: 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[],
893: PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX)
894: {
895: PetscDS prob;
896: PetscSF sf;
897: DM dmFace, dmCell, dmGrad;
898: const PetscScalar *facegeom, *cellgeom = NULL, *grad;
899: const PetscInt *leaves;
900: PetscScalar *x, *fx;
901: PetscInt dim, nleaves, loc, fStart, fEnd, pdim, i;
902: PetscErrorCode ierr, ierru = 0;
905: DMGetPointSF(dm, &sf);
906: PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
907: nleaves = PetscMax(0, nleaves);
908: DMGetDimension(dm, &dim);
909: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
910: DMGetDS(dm, &prob);
911: VecGetDM(faceGeometry, &dmFace);
912: VecGetArrayRead(faceGeometry, &facegeom);
913: if (cellGeometry) {
914: VecGetDM(cellGeometry, &dmCell);
915: VecGetArrayRead(cellGeometry, &cellgeom);
916: }
917: if (Grad) {
918: PetscFV fv;
920: PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);
921: VecGetDM(Grad, &dmGrad);
922: VecGetArrayRead(Grad, &grad);
923: PetscFVGetNumComponents(fv, &pdim);
924: DMGetWorkArray(dm, pdim, MPIU_SCALAR, &fx);
925: }
926: VecGetArray(locX, &x);
927: for (i = 0; i < numids; ++i) {
928: IS faceIS;
929: const PetscInt *faces;
930: PetscInt numFaces, f;
932: DMLabelGetStratumIS(label, ids[i], &faceIS);
933: if (!faceIS) continue; /* No points with that id on this process */
934: ISGetLocalSize(faceIS, &numFaces);
935: ISGetIndices(faceIS, &faces);
936: for (f = 0; f < numFaces; ++f) {
937: const PetscInt face = faces[f], *cells;
938: PetscFVFaceGeom *fg;
940: if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
941: PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);
942: if (loc >= 0) continue;
943: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
944: DMPlexGetSupport(dm, face, &cells);
945: if (Grad) {
946: PetscFVCellGeom *cg;
947: PetscScalar *cx, *cgrad;
948: PetscScalar *xG;
949: PetscReal dx[3];
950: PetscInt d;
952: DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);
953: DMPlexPointLocalRead(dm, cells[0], x, &cx);
954: DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);
955: DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
956: DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx);
957: for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx);
958: ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);
959: if (ierru) {
960: ISRestoreIndices(faceIS, &faces);
961: ISDestroy(&faceIS);
962: goto cleanup;
963: }
964: } else {
965: PetscScalar *xI;
966: PetscScalar *xG;
968: DMPlexPointLocalRead(dm, cells[0], x, &xI);
969: DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
970: ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);
971: if (ierru) {
972: ISRestoreIndices(faceIS, &faces);
973: ISDestroy(&faceIS);
974: goto cleanup;
975: }
976: }
977: }
978: ISRestoreIndices(faceIS, &faces);
979: ISDestroy(&faceIS);
980: }
981: cleanup:
982: VecRestoreArray(locX, &x);
983: if (Grad) {
984: DMRestoreWorkArray(dm, pdim, MPIU_SCALAR, &fx);
985: VecRestoreArrayRead(Grad, &grad);
986: }
987: if (cellGeometry) {VecRestoreArrayRead(cellGeometry, &cellgeom);}
988: VecRestoreArrayRead(faceGeometry, &facegeom);
989: CHKERRQ(ierru);
990: return(0);
991: }
993: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
994: {
995: PetscInt c;
996: for (c = 0; c < Nc; ++c) u[c] = 0.0;
997: return 0;
998: }
1000: PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
1001: {
1002: PetscObject isZero;
1003: PetscDS prob;
1004: PetscInt numBd, b;
1008: DMGetDS(dm, &prob);
1009: PetscDSGetNumBoundary(prob, &numBd);
1010: PetscObjectQuery((PetscObject) locX, "__Vec_bc_zero__", &isZero);
1011: for (b = 0; b < numBd; ++b) {
1012: DMBoundaryConditionType type;
1013: const char *name, *labelname;
1014: DMLabel label;
1015: PetscInt field, Nc;
1016: const PetscInt *comps;
1017: PetscObject obj;
1018: PetscClassId id;
1019: void (*func)(void);
1020: PetscInt numids;
1021: const PetscInt *ids;
1022: void *ctx;
1024: DMGetBoundary(dm, b, &type, &name, &labelname, &field, &Nc, &comps, &func, NULL, &numids, &ids, &ctx);
1025: if (insertEssential != (type & DM_BC_ESSENTIAL)) continue;
1026: DMGetLabel(dm, labelname, &label);
1027: if (!label) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "Label %s for boundary condition %s does not exist in the DM", labelname, name);
1028: DMGetField(dm, field, NULL, &obj);
1029: PetscObjectGetClassId(obj, &id);
1030: if (id == PETSCFE_CLASSID) {
1031: switch (type) {
1032: /* for FEM, there is no insertion to be done for non-essential boundary conditions */
1033: case DM_BC_ESSENTIAL:
1034: if (isZero) func = (void (*)(void)) zero;
1035: DMPlexLabelAddCells(dm,label);
1036: DMPlexInsertBoundaryValuesEssential(dm, time, field, Nc, comps, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);
1037: DMPlexLabelClearCells(dm,label);
1038: break;
1039: case DM_BC_ESSENTIAL_FIELD:
1040: DMPlexLabelAddCells(dm,label);
1041: DMPlexInsertBoundaryValuesEssentialField(dm, time, locX, field, Nc, comps, label, numids, ids,
1042: (void (*)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1043: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1044: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) func, ctx, locX);
1045: DMPlexLabelClearCells(dm,label);
1046: break;
1047: default: break;
1048: }
1049: } else if (id == PETSCFV_CLASSID) {
1050: if (!faceGeomFVM) continue;
1051: DMPlexInsertBoundaryValuesRiemann(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, field, Nc, comps, label, numids, ids,
1052: (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);
1053: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1054: }
1055: return(0);
1056: }
1058: PetscErrorCode DMPlexInsertTimeDerivativeBoundaryValues_Plex(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
1059: {
1060: PetscObject isZero;
1061: PetscDS prob;
1062: PetscInt numBd, b;
1066: if (!locX) return(0);
1067: DMGetDS(dm, &prob);
1068: PetscDSGetNumBoundary(prob, &numBd);
1069: PetscObjectQuery((PetscObject) locX, "__Vec_bc_zero__", &isZero);
1070: for (b = 0; b < numBd; ++b) {
1071: DMBoundaryConditionType type;
1072: const char *name, *labelname;
1073: DMLabel label;
1074: PetscInt field, Nc;
1075: const PetscInt *comps;
1076: PetscObject obj;
1077: PetscClassId id;
1078: void (*func_t)(void);
1079: PetscInt numids;
1080: const PetscInt *ids;
1081: void *ctx;
1083: DMGetBoundary(dm, b, &type, &name, &labelname, &field, &Nc, &comps, NULL, &func_t, &numids, &ids, &ctx);
1084: if (!func_t) continue;
1085: if (insertEssential != (type & DM_BC_ESSENTIAL)) continue;
1086: DMGetLabel(dm, labelname, &label);
1087: if (!label) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "Label %s for boundary condition %s does not exist in the DM", labelname, name);
1088: DMGetField(dm, field, NULL, &obj);
1089: PetscObjectGetClassId(obj, &id);
1090: if (id == PETSCFE_CLASSID) {
1091: switch (type) {
1092: /* for FEM, there is no insertion to be done for non-essential boundary conditions */
1093: case DM_BC_ESSENTIAL:
1094: if (isZero) func_t = (void (*)(void)) zero;
1095: DMPlexLabelAddCells(dm,label);
1096: DMPlexInsertBoundaryValuesEssential(dm, time, field, Nc, comps, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func_t, ctx, locX);
1097: DMPlexLabelClearCells(dm,label);
1098: break;
1099: case DM_BC_ESSENTIAL_FIELD:
1100: DMPlexLabelAddCells(dm,label);
1101: DMPlexInsertBoundaryValuesEssentialField(dm, time, locX, field, Nc, comps, label, numids, ids,
1102: (void (*)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1103: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1104: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) func_t, ctx, locX);
1105: DMPlexLabelClearCells(dm,label);
1106: break;
1107: default: break;
1108: }
1109: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1110: }
1111: return(0);
1112: }
1114: /*@
1115: DMPlexInsertBoundaryValues - Puts coefficients which represent boundary values into the local solution vector
1117: Input Parameters:
1118: + dm - The DM
1119: . insertEssential - Should I insert essential (e.g. Dirichlet) or inessential (e.g. Neumann) boundary conditions
1120: . time - The time
1121: . faceGeomFVM - Face geometry data for FV discretizations
1122: . cellGeomFVM - Cell geometry data for FV discretizations
1123: - gradFVM - Gradient reconstruction data for FV discretizations
1125: Output Parameters:
1126: . locX - Solution updated with boundary values
1128: Level: developer
1130: .seealso: DMProjectFunctionLabelLocal()
1131: @*/
1132: PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
1133: {
1142: PetscTryMethod(dm,"DMPlexInsertBoundaryValues_C",(DM,PetscBool,Vec,PetscReal,Vec,Vec,Vec),(dm,insertEssential,locX,time,faceGeomFVM,cellGeomFVM,gradFVM));
1143: return(0);
1144: }
1146: /*@
1147: DMPlexInsertTimeDerivativeBoundaryValues - Puts coefficients which represent boundary values of the time derviative into the local solution vector
1149: Input Parameters:
1150: + dm - The DM
1151: . insertEssential - Should I insert essential (e.g. Dirichlet) or inessential (e.g. Neumann) boundary conditions
1152: . time - The time
1153: . faceGeomFVM - Face geometry data for FV discretizations
1154: . cellGeomFVM - Cell geometry data for FV discretizations
1155: - gradFVM - Gradient reconstruction data for FV discretizations
1157: Output Parameters:
1158: . locX_t - Solution updated with boundary values
1160: Level: developer
1162: .seealso: DMProjectFunctionLabelLocal()
1163: @*/
1164: PetscErrorCode DMPlexInsertTimeDerivativeBoundaryValues(DM dm, PetscBool insertEssential, Vec locX_t, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
1165: {
1174: PetscTryMethod(dm,"DMPlexInsertTimeDerviativeBoundaryValues_C",(DM,PetscBool,Vec,PetscReal,Vec,Vec,Vec),(dm,insertEssential,locX_t,time,faceGeomFVM,cellGeomFVM,gradFVM));
1175: return(0);
1176: }
1178: PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
1179: {
1180: Vec localX;
1181: PetscErrorCode ierr;
1184: DMGetLocalVector(dm, &localX);
1185: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, time, NULL, NULL, NULL);
1186: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1187: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1188: DMPlexComputeL2DiffLocal(dm, time, funcs, ctxs, localX, diff);
1189: DMRestoreLocalVector(dm, &localX);
1190: return(0);
1191: }
1193: /*@C
1194: DMComputeL2DiffLocal - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
1196: Collective on dm
1198: Input Parameters:
1199: + dm - The DM
1200: . time - The time
1201: . funcs - The functions to evaluate for each field component
1202: . ctxs - Optional array of contexts to pass to each function, or NULL.
1203: - localX - The coefficient vector u_h, a local vector
1205: Output Parameter:
1206: . diff - The diff ||u - u_h||_2
1208: Level: developer
1210: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
1211: @*/
1212: PetscErrorCode DMPlexComputeL2DiffLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec localX, PetscReal *diff)
1213: {
1214: const PetscInt debug = ((DM_Plex*)dm->data)->printL2;
1215: DM tdm;
1216: Vec tv;
1217: PetscSection section;
1218: PetscQuadrature quad;
1219: PetscFEGeom fegeom;
1220: PetscScalar *funcVal, *interpolant;
1221: PetscReal *coords, *gcoords;
1222: PetscReal localDiff = 0.0;
1223: const PetscReal *quadWeights;
1224: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cellHeight, cStart, cEnd, c, field, fieldOffset;
1225: PetscBool transform;
1226: PetscErrorCode ierr;
1229: DMGetDimension(dm, &dim);
1230: DMGetCoordinateDim(dm, &coordDim);
1231: fegeom.dimEmbed = coordDim;
1232: DMGetLocalSection(dm, §ion);
1233: PetscSectionGetNumFields(section, &numFields);
1234: DMGetBasisTransformDM_Internal(dm, &tdm);
1235: DMGetBasisTransformVec_Internal(dm, &tv);
1236: DMHasBasisTransform(dm, &transform);
1237: for (field = 0; field < numFields; ++field) {
1238: PetscObject obj;
1239: PetscClassId id;
1240: PetscInt Nc;
1242: DMGetField(dm, field, NULL, &obj);
1243: PetscObjectGetClassId(obj, &id);
1244: if (id == PETSCFE_CLASSID) {
1245: PetscFE fe = (PetscFE) obj;
1247: PetscFEGetQuadrature(fe, &quad);
1248: PetscFEGetNumComponents(fe, &Nc);
1249: } else if (id == PETSCFV_CLASSID) {
1250: PetscFV fv = (PetscFV) obj;
1252: PetscFVGetQuadrature(fv, &quad);
1253: PetscFVGetNumComponents(fv, &Nc);
1254: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1255: numComponents += Nc;
1256: }
1257: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
1258: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1259: PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);
1260: DMPlexGetVTKCellHeight(dm, &cellHeight);
1261: DMPlexGetSimplexOrBoxCells(dm, cellHeight, &cStart, &cEnd);
1262: for (c = cStart; c < cEnd; ++c) {
1263: PetscScalar *x = NULL;
1264: PetscReal elemDiff = 0.0;
1265: PetscInt qc = 0;
1267: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
1268: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
1270: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1271: PetscObject obj;
1272: PetscClassId id;
1273: void * const ctx = ctxs ? ctxs[field] : NULL;
1274: PetscInt Nb, Nc, q, fc;
1276: DMGetField(dm, field, NULL, &obj);
1277: PetscObjectGetClassId(obj, &id);
1278: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1279: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1280: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1281: if (debug) {
1282: char title[1024];
1283: PetscSNPrintf(title, 1023, "Solution for Field %D", field);
1284: DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
1285: }
1286: for (q = 0; q < Nq; ++q) {
1287: PetscFEGeom qgeom;
1289: qgeom.dimEmbed = fegeom.dimEmbed;
1290: qgeom.J = &fegeom.J[q*coordDim*coordDim];
1291: qgeom.invJ = &fegeom.invJ[q*coordDim*coordDim];
1292: qgeom.detJ = &fegeom.detJ[q];
1293: 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);
1294: if (transform) {
1295: gcoords = &coords[coordDim*Nq];
1296: DMPlexBasisTransformApplyReal_Internal(dm, &coords[coordDim*q], PETSC_TRUE, coordDim, &coords[coordDim*q], gcoords, dm->transformCtx);
1297: } else {
1298: gcoords = &coords[coordDim*q];
1299: }
1300: (*funcs[field])(coordDim, time, gcoords, Nc, funcVal, ctx);
1301: if (ierr) {
1302: PetscErrorCode ierr2;
1303: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
1304: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1305: ierr2 = PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2);
1306:
1307: }
1308: if (transform) {DMPlexBasisTransformApply_Internal(dm, &coords[coordDim*q], PETSC_FALSE, Nc, funcVal, funcVal, dm->transformCtx);}
1309: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], &qgeom, q, interpolant);}
1310: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
1311: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1312: for (fc = 0; fc < Nc; ++fc) {
1313: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
1314: if (debug) {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]));}
1315: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q];
1316: }
1317: }
1318: fieldOffset += Nb;
1319: qc += Nc;
1320: }
1321: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1322: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %D diff %g\n", c, (double)elemDiff);}
1323: localDiff += elemDiff;
1324: }
1325: PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
1326: MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1327: *diff = PetscSqrtReal(*diff);
1328: return(0);
1329: }
1331: 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)
1332: {
1333: const PetscInt debug = ((DM_Plex*)dm->data)->printL2;
1334: DM tdm;
1335: PetscSection section;
1336: PetscQuadrature quad;
1337: Vec localX, tv;
1338: PetscScalar *funcVal, *interpolant;
1339: const PetscReal *quadWeights;
1340: PetscFEGeom fegeom;
1341: PetscReal *coords, *gcoords;
1342: PetscReal localDiff = 0.0;
1343: PetscInt dim, coordDim, qNc = 0, Nq = 0, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset;
1344: PetscBool transform;
1345: PetscErrorCode ierr;
1348: DMGetDimension(dm, &dim);
1349: DMGetCoordinateDim(dm, &coordDim);
1350: fegeom.dimEmbed = coordDim;
1351: DMGetLocalSection(dm, §ion);
1352: PetscSectionGetNumFields(section, &numFields);
1353: DMGetLocalVector(dm, &localX);
1354: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1355: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1356: DMGetBasisTransformDM_Internal(dm, &tdm);
1357: DMGetBasisTransformVec_Internal(dm, &tv);
1358: DMHasBasisTransform(dm, &transform);
1359: for (field = 0; field < numFields; ++field) {
1360: PetscFE fe;
1361: PetscInt Nc;
1363: DMGetField(dm, field, NULL, (PetscObject *) &fe);
1364: PetscFEGetQuadrature(fe, &quad);
1365: PetscFEGetNumComponents(fe, &Nc);
1366: numComponents += Nc;
1367: }
1368: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
1369: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1370: /* DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
1371: PetscMalloc6(numComponents,&funcVal,coordDim*Nq,&coords,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ,numComponents*coordDim,&interpolant,Nq,&fegeom.detJ);
1372: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
1373: for (c = cStart; c < cEnd; ++c) {
1374: PetscScalar *x = NULL;
1375: PetscReal elemDiff = 0.0;
1376: PetscInt qc = 0;
1378: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
1379: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
1381: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1382: PetscFE fe;
1383: void * const ctx = ctxs ? ctxs[field] : NULL;
1384: PetscInt Nb, Nc, q, fc;
1386: DMGetField(dm, field, NULL, (PetscObject *) &fe);
1387: PetscFEGetDimension(fe, &Nb);
1388: PetscFEGetNumComponents(fe, &Nc);
1389: if (debug) {
1390: char title[1024];
1391: PetscSNPrintf(title, 1023, "Solution for Field %D", field);
1392: DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
1393: }
1394: for (q = 0; q < Nq; ++q) {
1395: PetscFEGeom qgeom;
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 points %D", (double)fegeom.detJ[q], c, q);
1402: if (transform) {
1403: gcoords = &coords[coordDim*Nq];
1404: DMPlexBasisTransformApplyReal_Internal(dm, &coords[coordDim*q], PETSC_TRUE, coordDim, &coords[coordDim*q], gcoords, dm->transformCtx);
1405: } else {
1406: gcoords = &coords[coordDim*q];
1407: }
1408: (*funcs[field])(coordDim, time, gcoords, n, 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 = PetscFree6(funcVal,coords,fegeom.J,fegeom.invJ,interpolant,fegeom.detJ);CHKERRQ(ierr2);
1414:
1415: }
1416: if (transform) {DMPlexBasisTransformApply_Internal(dm, &coords[coordDim*q], PETSC_FALSE, Nc, funcVal, funcVal, dm->transformCtx);}
1417: PetscFEInterpolateGradient_Static(fe, &x[fieldOffset], &qgeom, q, interpolant);
1418: /* Overwrite with the dot product if the normal is given */
1419: if (n) {
1420: for (fc = 0; fc < Nc; ++fc) {
1421: PetscScalar sum = 0.0;
1422: PetscInt d;
1423: for (d = 0; d < dim; ++d) sum += interpolant[fc*dim+d]*n[d];
1424: interpolant[fc] = sum;
1425: }
1426: }
1427: for (fc = 0; fc < Nc; ++fc) {
1428: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
1429: if (debug) {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]));}
1430: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q];
1431: }
1432: }
1433: fieldOffset += Nb;
1434: qc += Nc;
1435: }
1436: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1437: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %D diff %g\n", c, (double)elemDiff);}
1438: localDiff += elemDiff;
1439: }
1440: PetscFree6(funcVal,coords,fegeom.J,fegeom.invJ,interpolant,fegeom.detJ);
1441: DMRestoreLocalVector(dm, &localX);
1442: MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1443: *diff = PetscSqrtReal(*diff);
1444: return(0);
1445: }
1447: PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
1448: {
1449: const PetscInt debug = ((DM_Plex*)dm->data)->printL2;
1450: DM tdm;
1451: DMLabel depthLabel;
1452: PetscSection section;
1453: Vec localX, tv;
1454: PetscReal *localDiff;
1455: PetscInt dim, depth, dE, Nf, f, Nds, s;
1456: PetscBool transform;
1457: PetscErrorCode ierr;
1460: DMGetDimension(dm, &dim);
1461: DMGetCoordinateDim(dm, &dE);
1462: DMGetLocalSection(dm, §ion);
1463: DMGetLocalVector(dm, &localX);
1464: DMGetBasisTransformDM_Internal(dm, &tdm);
1465: DMGetBasisTransformVec_Internal(dm, &tv);
1466: DMHasBasisTransform(dm, &transform);
1467: DMGetNumFields(dm, &Nf);
1468: DMPlexGetDepthLabel(dm, &depthLabel);
1469: DMLabelGetNumValues(depthLabel, &depth);
1471: VecSet(localX, 0.0);
1472: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1473: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1474: DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
1475: DMGetNumDS(dm, &Nds);
1476: PetscCalloc1(Nf, &localDiff);
1477: for (s = 0; s < Nds; ++s) {
1478: PetscDS ds;
1479: DMLabel label;
1480: IS fieldIS, pointIS;
1481: const PetscInt *fields, *points = NULL;
1482: PetscQuadrature quad;
1483: const PetscReal *quadPoints, *quadWeights;
1484: PetscFEGeom fegeom;
1485: PetscReal *coords, *gcoords;
1486: PetscScalar *funcVal, *interpolant;
1487: PetscBool isHybrid;
1488: PetscInt qNc, Nq, totNc, cStart = 0, cEnd, c, dsNf;
1490: DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds);
1491: ISGetIndices(fieldIS, &fields);
1492: PetscDSGetHybrid(ds, &isHybrid);
1493: PetscDSGetNumFields(ds, &dsNf);
1494: PetscDSGetTotalComponents(ds, &totNc);
1495: PetscDSGetQuadrature(ds, &quad);
1496: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1497: if ((qNc != 1) && (qNc != totNc)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, totNc);
1498: PetscCalloc6(totNc, &funcVal, totNc, &interpolant, dE*(Nq+1), &coords,Nq, &fegeom.detJ, dE*dE*Nq, &fegeom.J, dE*dE*Nq, &fegeom.invJ);
1499: if (!label) {
1500: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
1501: } else {
1502: DMLabelGetStratumIS(label, 1, &pointIS);
1503: ISGetLocalSize(pointIS, &cEnd);
1504: ISGetIndices(pointIS, &points);
1505: }
1506: for (c = cStart; c < cEnd; ++c) {
1507: const PetscInt cell = points ? points[c] : c;
1508: PetscScalar *x = NULL;
1509: PetscInt qc = 0, fOff = 0, dep, fStart = isHybrid ? dsNf-1 : 0;
1511: DMLabelGetValue(depthLabel, cell, &dep);
1512: if (dep != depth-1) continue;
1513: if (isHybrid) {
1514: const PetscInt *cone;
1516: DMPlexGetCone(dm, cell, &cone);
1517: DMPlexComputeCellGeometryFEM(dm, cone[0], quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
1518: } else {
1519: DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
1520: }
1521: DMPlexVecGetClosure(dm, NULL, localX, cell, NULL, &x);
1522: for (f = fStart; f < dsNf; ++f) {
1523: PetscObject obj;
1524: PetscClassId id;
1525: void * const ctx = ctxs ? ctxs[fields[f]] : NULL;
1526: PetscInt Nb, Nc, q, fc;
1527: PetscReal elemDiff = 0.0;
1529: PetscDSGetDiscretization(ds, f, &obj);
1530: PetscObjectGetClassId(obj, &id);
1531: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1532: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1533: else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", fields[f]);
1534: if (debug) {
1535: char title[1024];
1536: PetscSNPrintf(title, 1023, "Solution for Field %D", fields[f]);
1537: DMPrintCellVector(cell, title, Nb, &x[fOff]);
1538: }
1539: for (q = 0; q < Nq; ++q) {
1540: PetscFEGeom qgeom;
1542: qgeom.dimEmbed = fegeom.dimEmbed;
1543: qgeom.J = &fegeom.J[q*dE*dE];
1544: qgeom.invJ = &fegeom.invJ[q*dE*dE];
1545: qgeom.detJ = &fegeom.detJ[q];
1546: if (fegeom.detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for cell %D, quadrature point %D", (double)fegeom.detJ[q], cell, q);
1547: if (transform) {
1548: gcoords = &coords[dE*Nq];
1549: DMPlexBasisTransformApplyReal_Internal(dm, &coords[dE*q], PETSC_TRUE, dE, &coords[dE*q], gcoords, dm->transformCtx);
1550: } else {
1551: gcoords = &coords[dE*q];
1552: }
1553: (*funcs[fields[f]])(dE, time, gcoords, Nc, funcVal, ctx);
1554: if (ierr) {
1555: PetscErrorCode ierr2;
1556: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, cell, NULL, &x);CHKERRQ(ierr2);
1557: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1558: ierr2 = PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2);
1559:
1560: }
1561: if (transform) {DMPlexBasisTransformApply_Internal(dm, &coords[dE*q], PETSC_FALSE, Nc, funcVal, funcVal, dm->transformCtx);}
1562: /* Call once for each face, except for lagrange field */
1563: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fOff], &qgeom, q, interpolant);}
1564: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fOff], q, interpolant);}
1565: else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", fields[f]);
1566: for (fc = 0; fc < Nc; ++fc) {
1567: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
1568: if (debug) {PetscPrintf(PETSC_COMM_SELF, " cell %D field %D,%D point %g %g %g diff %g\n", cell, fields[f], fc, (double)(dE > 0 ? coords[dE*q] : 0.), (double)(dE > 1 ? coords[dE*q+1] : 0.),(double)(dE > 2 ? coords[dE*q+2] : 0.), (double)(PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q]));}
1569: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q];
1570: }
1571: }
1572: fOff += Nb;
1573: qc += Nc;
1574: localDiff[fields[f]] += elemDiff;
1575: if (debug) {PetscPrintf(PETSC_COMM_SELF, " cell %D field %D cum diff %g\n", cell, fields[f], (double)localDiff[fields[f]]);}
1576: }
1577: DMPlexVecRestoreClosure(dm, NULL, localX, cell, NULL, &x);
1578: }
1579: if (label) {
1580: ISRestoreIndices(pointIS, &points);
1581: ISDestroy(&pointIS);
1582: }
1583: ISRestoreIndices(fieldIS, &fields);
1584: PetscFree6(funcVal, interpolant, coords, fegeom.detJ, fegeom.J, fegeom.invJ);
1585: }
1586: DMRestoreLocalVector(dm, &localX);
1587: MPIU_Allreduce(localDiff, diff, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1588: PetscFree(localDiff);
1589: for (f = 0; f < Nf; ++f) diff[f] = PetscSqrtReal(diff[f]);
1590: return(0);
1591: }
1593: /*@C
1594: 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.
1596: Collective on dm
1598: Input Parameters:
1599: + dm - The DM
1600: . time - The time
1601: . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation
1602: . ctxs - Optional array of contexts to pass to each function, or NULL.
1603: - X - The coefficient vector u_h
1605: Output Parameter:
1606: . D - A Vec which holds the difference ||u - u_h||_2 for each cell
1608: Level: developer
1610: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
1611: @*/
1612: PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D)
1613: {
1614: PetscSection section;
1615: PetscQuadrature quad;
1616: Vec localX;
1617: PetscFEGeom fegeom;
1618: PetscScalar *funcVal, *interpolant;
1619: PetscReal *coords;
1620: const PetscReal *quadPoints, *quadWeights;
1621: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, c, field, fieldOffset;
1622: PetscErrorCode ierr;
1625: VecSet(D, 0.0);
1626: DMGetDimension(dm, &dim);
1627: DMGetCoordinateDim(dm, &coordDim);
1628: DMGetLocalSection(dm, §ion);
1629: PetscSectionGetNumFields(section, &numFields);
1630: DMGetLocalVector(dm, &localX);
1631: DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
1632: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1633: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1634: for (field = 0; field < numFields; ++field) {
1635: PetscObject obj;
1636: PetscClassId id;
1637: PetscInt Nc;
1639: DMGetField(dm, field, NULL, &obj);
1640: PetscObjectGetClassId(obj, &id);
1641: if (id == PETSCFE_CLASSID) {
1642: PetscFE fe = (PetscFE) obj;
1644: PetscFEGetQuadrature(fe, &quad);
1645: PetscFEGetNumComponents(fe, &Nc);
1646: } else if (id == PETSCFV_CLASSID) {
1647: PetscFV fv = (PetscFV) obj;
1649: PetscFVGetQuadrature(fv, &quad);
1650: PetscFVGetNumComponents(fv, &Nc);
1651: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1652: numComponents += Nc;
1653: }
1654: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1655: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1656: PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);
1657: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
1658: for (c = cStart; c < cEnd; ++c) {
1659: PetscScalar *x = NULL;
1660: PetscScalar elemDiff = 0.0;
1661: PetscInt qc = 0;
1663: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
1664: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
1666: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1667: PetscObject obj;
1668: PetscClassId id;
1669: void * const ctx = ctxs ? ctxs[field] : NULL;
1670: PetscInt Nb, Nc, q, fc;
1672: DMGetField(dm, field, NULL, &obj);
1673: PetscObjectGetClassId(obj, &id);
1674: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1675: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1676: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1677: if (funcs[field]) {
1678: for (q = 0; q < Nq; ++q) {
1679: PetscFEGeom qgeom;
1681: qgeom.dimEmbed = fegeom.dimEmbed;
1682: qgeom.J = &fegeom.J[q*coordDim*coordDim];
1683: qgeom.invJ = &fegeom.invJ[q*coordDim*coordDim];
1684: qgeom.detJ = &fegeom.detJ[q];
1685: 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);
1686: (*funcs[field])(coordDim, time, &coords[q*coordDim], Nc, funcVal, ctx);
1687: if (ierr) {
1688: PetscErrorCode ierr2;
1689: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
1690: ierr2 = PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2);
1691: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1692:
1693: }
1694: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], &qgeom, q, interpolant);}
1695: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
1696: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1697: for (fc = 0; fc < Nc; ++fc) {
1698: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
1699: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q];
1700: }
1701: }
1702: }
1703: fieldOffset += Nb;
1704: qc += Nc;
1705: }
1706: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1707: VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);
1708: }
1709: PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
1710: DMRestoreLocalVector(dm, &localX);
1711: VecSqrtAbs(D);
1712: return(0);
1713: }
1715: /*@C
1716: DMPlexComputeGradientClementInterpolant - This function computes the L2 projection of the cellwise gradient of a function u onto P1, and stores it in a Vec.
1718: Collective on dm
1720: Input Parameters:
1721: + dm - The DM
1722: - LocX - The coefficient vector u_h
1724: Output Parameter:
1725: . locC - A Vec which holds the Clement interpolant of the gradient
1727: Notes:
1728: Add citation to (Clement, 1975) and definition of the interpolant
1729: \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
1731: Level: developer
1733: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
1734: @*/
1735: PetscErrorCode DMPlexComputeGradientClementInterpolant(DM dm, Vec locX, Vec locC)
1736: {
1737: DM_Plex *mesh = (DM_Plex *) dm->data;
1738: PetscInt debug = mesh->printFEM;
1739: DM dmC;
1740: PetscSection section;
1741: PetscQuadrature quad;
1742: PetscScalar *interpolant, *gradsum;
1743: PetscFEGeom fegeom;
1744: PetscReal *coords;
1745: const PetscReal *quadPoints, *quadWeights;
1746: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, vStart, vEnd, v, field, fieldOffset;
1747: PetscErrorCode ierr;
1750: VecGetDM(locC, &dmC);
1751: VecSet(locC, 0.0);
1752: DMGetDimension(dm, &dim);
1753: DMGetCoordinateDim(dm, &coordDim);
1754: fegeom.dimEmbed = coordDim;
1755: DMGetLocalSection(dm, §ion);
1756: PetscSectionGetNumFields(section, &numFields);
1757: for (field = 0; field < numFields; ++field) {
1758: PetscObject obj;
1759: PetscClassId id;
1760: PetscInt Nc;
1762: DMGetField(dm, field, NULL, &obj);
1763: PetscObjectGetClassId(obj, &id);
1764: if (id == PETSCFE_CLASSID) {
1765: PetscFE fe = (PetscFE) obj;
1767: PetscFEGetQuadrature(fe, &quad);
1768: PetscFEGetNumComponents(fe, &Nc);
1769: } else if (id == PETSCFV_CLASSID) {
1770: PetscFV fv = (PetscFV) obj;
1772: PetscFVGetQuadrature(fv, &quad);
1773: PetscFVGetNumComponents(fv, &Nc);
1774: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1775: numComponents += Nc;
1776: }
1777: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1778: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1779: PetscMalloc6(coordDim*numComponents*2,&gradsum,coordDim*numComponents,&interpolant,coordDim*Nq,&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);
1780: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1781: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
1782: for (v = vStart; v < vEnd; ++v) {
1783: PetscScalar volsum = 0.0;
1784: PetscInt *star = NULL;
1785: PetscInt starSize, st, d, fc;
1787: PetscArrayzero(gradsum, coordDim*numComponents);
1788: DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1789: for (st = 0; st < starSize*2; st += 2) {
1790: const PetscInt cell = star[st];
1791: PetscScalar *grad = &gradsum[coordDim*numComponents];
1792: PetscScalar *x = NULL;
1793: PetscReal vol = 0.0;
1795: if ((cell < cStart) || (cell >= cEnd)) continue;
1796: DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
1797: DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x);
1798: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1799: PetscObject obj;
1800: PetscClassId id;
1801: PetscInt Nb, Nc, q, qc = 0;
1803: PetscArrayzero(grad, coordDim*numComponents);
1804: DMGetField(dm, field, NULL, &obj);
1805: PetscObjectGetClassId(obj, &id);
1806: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1807: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1808: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1809: for (q = 0; q < Nq; ++q) {
1810: PetscFEGeom qgeom;
1812: qgeom.dimEmbed = fegeom.dimEmbed;
1813: qgeom.J = &fegeom.J[q*coordDim*coordDim];
1814: qgeom.invJ = &fegeom.invJ[q*coordDim*coordDim];
1815: qgeom.detJ = &fegeom.detJ[q];
1816: 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);
1817: if (ierr) {
1818: PetscErrorCode ierr2;
1819: ierr2 = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr2);
1820: ierr2 = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr2);
1821: ierr2 = PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2);
1822:
1823: }
1824: if (id == PETSCFE_CLASSID) {PetscFEInterpolateGradient_Static((PetscFE) obj, &x[fieldOffset], &qgeom, q, interpolant);}
1825: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1826: for (fc = 0; fc < Nc; ++fc) {
1827: const PetscReal wt = quadWeights[q*qNc+qc+fc];
1829: for (d = 0; d < coordDim; ++d) grad[fc*coordDim+d] += interpolant[fc*dim+d]*wt*fegeom.detJ[q];
1830: }
1831: vol += quadWeights[q*qNc]*fegeom.detJ[q];
1832: }
1833: fieldOffset += Nb;
1834: qc += Nc;
1835: }
1836: DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);
1837: for (fc = 0; fc < numComponents; ++fc) {
1838: for (d = 0; d < coordDim; ++d) {
1839: gradsum[fc*coordDim+d] += grad[fc*coordDim+d];
1840: }
1841: }
1842: volsum += vol;
1843: if (debug) {
1844: PetscPrintf(PETSC_COMM_SELF, "Cell %D gradient: [", cell);
1845: for (fc = 0; fc < numComponents; ++fc) {
1846: for (d = 0; d < coordDim; ++d) {
1847: if (fc || d > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
1848: PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc*coordDim+d]));
1849: }
1850: }
1851: PetscPrintf(PETSC_COMM_SELF, "]\n");
1852: }
1853: }
1854: for (fc = 0; fc < numComponents; ++fc) {
1855: for (d = 0; d < coordDim; ++d) gradsum[fc*coordDim+d] /= volsum;
1856: }
1857: DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1858: DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES);
1859: }
1860: PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
1861: return(0);
1862: }
1864: static PetscErrorCode DMPlexComputeIntegral_Internal(DM dm, Vec X, PetscInt cStart, PetscInt cEnd, PetscScalar *cintegral, void *user)
1865: {
1866: DM dmAux = NULL;
1867: PetscDS prob, probAux = NULL;
1868: PetscSection section, sectionAux;
1869: Vec locX, locA;
1870: PetscInt dim, numCells = cEnd - cStart, c, f;
1871: PetscBool useFVM = PETSC_FALSE;
1872: /* DS */
1873: PetscInt Nf, totDim, *uOff, *uOff_x, numConstants;
1874: PetscInt NfAux, totDimAux, *aOff;
1875: PetscScalar *u, *a;
1876: const PetscScalar *constants;
1877: /* Geometry */
1878: PetscFEGeom *cgeomFEM;
1879: DM dmGrad;
1880: PetscQuadrature affineQuad = NULL;
1881: Vec cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL;
1882: PetscFVCellGeom *cgeomFVM;
1883: const PetscScalar *lgrad;
1884: PetscInt maxDegree;
1885: DMField coordField;
1886: IS cellIS;
1887: PetscErrorCode ierr;
1890: DMGetDS(dm, &prob);
1891: DMGetDimension(dm, &dim);
1892: DMGetLocalSection(dm, §ion);
1893: PetscSectionGetNumFields(section, &Nf);
1894: /* Determine which discretizations we have */
1895: for (f = 0; f < Nf; ++f) {
1896: PetscObject obj;
1897: PetscClassId id;
1899: PetscDSGetDiscretization(prob, f, &obj);
1900: PetscObjectGetClassId(obj, &id);
1901: if (id == PETSCFV_CLASSID) useFVM = PETSC_TRUE;
1902: }
1903: /* Get local solution with boundary values */
1904: DMGetLocalVector(dm, &locX);
1905: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);
1906: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
1907: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
1908: /* Read DS information */
1909: PetscDSGetTotalDimension(prob, &totDim);
1910: PetscDSGetComponentOffsets(prob, &uOff);
1911: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
1912: ISCreateStride(PETSC_COMM_SELF,numCells,cStart,1,&cellIS);
1913: PetscDSGetConstants(prob, &numConstants, &constants);
1914: /* Read Auxiliary DS information */
1915: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1916: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
1917: if (dmAux) {
1918: DMGetDS(dmAux, &probAux);
1919: PetscDSGetNumFields(probAux, &NfAux);
1920: DMGetLocalSection(dmAux, §ionAux);
1921: PetscDSGetTotalDimension(probAux, &totDimAux);
1922: PetscDSGetComponentOffsets(probAux, &aOff);
1923: }
1924: /* Allocate data arrays */
1925: PetscCalloc1(numCells*totDim, &u);
1926: if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1927: /* Read out geometry */
1928: DMGetCoordinateField(dm,&coordField);
1929: DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);
1930: if (maxDegree <= 1) {
1931: DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);
1932: if (affineQuad) {
1933: DMFieldCreateFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&cgeomFEM);
1934: }
1935: }
1936: if (useFVM) {
1937: PetscFV fv = NULL;
1938: Vec grad;
1939: PetscInt fStart, fEnd;
1940: PetscBool compGrad;
1942: for (f = 0; f < Nf; ++f) {
1943: PetscObject obj;
1944: PetscClassId id;
1946: PetscDSGetDiscretization(prob, f, &obj);
1947: PetscObjectGetClassId(obj, &id);
1948: if (id == PETSCFV_CLASSID) {fv = (PetscFV) obj; break;}
1949: }
1950: PetscFVGetComputeGradients(fv, &compGrad);
1951: PetscFVSetComputeGradients(fv, PETSC_TRUE);
1952: DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);
1953: DMPlexComputeGradientFVM(dm, fv, faceGeometryFVM, cellGeometryFVM, &dmGrad);
1954: PetscFVSetComputeGradients(fv, compGrad);
1955: VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1956: /* Reconstruct and limit cell gradients */
1957: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1958: DMGetGlobalVector(dmGrad, &grad);
1959: DMPlexReconstructGradients_Internal(dm, fv, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
1960: /* Communicate gradient values */
1961: DMGetLocalVector(dmGrad, &locGrad);
1962: DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);
1963: DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);
1964: DMRestoreGlobalVector(dmGrad, &grad);
1965: /* Handle non-essential (e.g. outflow) boundary values */
1966: DMPlexInsertBoundaryValues(dm, PETSC_FALSE, locX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);
1967: VecGetArrayRead(locGrad, &lgrad);
1968: }
1969: /* Read out data from inputs */
1970: for (c = cStart; c < cEnd; ++c) {
1971: PetscScalar *x = NULL;
1972: PetscInt i;
1974: DMPlexVecGetClosure(dm, section, locX, c, NULL, &x);
1975: for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1976: DMPlexVecRestoreClosure(dm, section, locX, c, NULL, &x);
1977: if (dmAux) {
1978: DMPlexVecGetClosure(dmAux, sectionAux, locA, c, NULL, &x);
1979: for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1980: DMPlexVecRestoreClosure(dmAux, sectionAux, locA, c, NULL, &x);
1981: }
1982: }
1983: /* Do integration for each field */
1984: for (f = 0; f < Nf; ++f) {
1985: PetscObject obj;
1986: PetscClassId id;
1987: PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;
1989: PetscDSGetDiscretization(prob, f, &obj);
1990: PetscObjectGetClassId(obj, &id);
1991: if (id == PETSCFE_CLASSID) {
1992: PetscFE fe = (PetscFE) obj;
1993: PetscQuadrature q;
1994: PetscFEGeom *chunkGeom = NULL;
1995: PetscInt Nq, Nb;
1997: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1998: PetscFEGetQuadrature(fe, &q);
1999: PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);
2000: PetscFEGetDimension(fe, &Nb);
2001: blockSize = Nb*Nq;
2002: batchSize = numBlocks * blockSize;
2003: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2004: numChunks = numCells / (numBatches*batchSize);
2005: Ne = numChunks*numBatches*batchSize;
2006: Nr = numCells % (numBatches*batchSize);
2007: offset = numCells - Nr;
2008: if (!affineQuad) {
2009: DMFieldCreateFEGeom(coordField,cellIS,q,PETSC_FALSE,&cgeomFEM);
2010: }
2011: PetscFEGeomGetChunk(cgeomFEM,0,offset,&chunkGeom);
2012: PetscFEIntegrate(prob, f, Ne, chunkGeom, u, probAux, a, cintegral);
2013: PetscFEGeomGetChunk(cgeomFEM,offset,numCells,&chunkGeom);
2014: PetscFEIntegrate(prob, f, Nr, chunkGeom, &u[offset*totDim], probAux, &a[offset*totDimAux], &cintegral[offset*Nf]);
2015: PetscFEGeomRestoreChunk(cgeomFEM,offset,numCells,&chunkGeom);
2016: if (!affineQuad) {
2017: PetscFEGeomDestroy(&cgeomFEM);
2018: }
2019: } else if (id == PETSCFV_CLASSID) {
2020: PetscInt foff;
2021: PetscPointFunc obj_func;
2022: PetscScalar lint;
2024: PetscDSGetObjective(prob, f, &obj_func);
2025: PetscDSGetFieldOffset(prob, f, &foff);
2026: if (obj_func) {
2027: for (c = 0; c < numCells; ++c) {
2028: PetscScalar *u_x;
2030: DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);
2031: 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);
2032: cintegral[c*Nf+f] += PetscRealPart(lint)*cgeomFVM[c].volume;
2033: }
2034: }
2035: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f);
2036: }
2037: /* Cleanup data arrays */
2038: if (useFVM) {
2039: VecRestoreArrayRead(locGrad, &lgrad);
2040: VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
2041: DMRestoreLocalVector(dmGrad, &locGrad);
2042: VecDestroy(&faceGeometryFVM);
2043: VecDestroy(&cellGeometryFVM);
2044: DMDestroy(&dmGrad);
2045: }
2046: if (dmAux) {PetscFree(a);}
2047: PetscFree(u);
2048: /* Cleanup */
2049: if (affineQuad) {
2050: PetscFEGeomDestroy(&cgeomFEM);
2051: }
2052: PetscQuadratureDestroy(&affineQuad);
2053: ISDestroy(&cellIS);
2054: DMRestoreLocalVector(dm, &locX);
2055: return(0);
2056: }
2058: /*@
2059: DMPlexComputeIntegralFEM - Form the integral over the domain from the global input X using pointwise functions specified by the user
2061: Input Parameters:
2062: + dm - The mesh
2063: . X - Global input vector
2064: - user - The user context
2066: Output Parameter:
2067: . integral - Integral for each field
2069: Level: developer
2071: .seealso: DMPlexComputeResidualFEM()
2072: @*/
2073: PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscScalar *integral, void *user)
2074: {
2075: DM_Plex *mesh = (DM_Plex *) dm->data;
2076: PetscScalar *cintegral, *lintegral;
2077: PetscInt Nf, f, cellHeight, cStart, cEnd, cell;
2084: PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
2085: DMGetNumFields(dm, &Nf);
2086: DMPlexGetVTKCellHeight(dm, &cellHeight);
2087: DMPlexGetSimplexOrBoxCells(dm, cellHeight, &cStart, &cEnd);
2088: /* TODO Introduce a loop over large chunks (right now this is a single chunk) */
2089: PetscCalloc2(Nf, &lintegral, (cEnd-cStart)*Nf, &cintegral);
2090: DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);
2091: /* Sum up values */
2092: for (cell = cStart; cell < cEnd; ++cell) {
2093: const PetscInt c = cell - cStart;
2095: if (mesh->printFEM > 1) {DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);}
2096: for (f = 0; f < Nf; ++f) lintegral[f] += cintegral[c*Nf+f];
2097: }
2098: MPIU_Allreduce(lintegral, integral, Nf, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject) dm));
2099: if (mesh->printFEM) {
2100: PetscPrintf(PetscObjectComm((PetscObject) dm), "Integral:");
2101: for (f = 0; f < Nf; ++f) {PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", (double) PetscRealPart(integral[f]));}
2102: PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");
2103: }
2104: PetscFree2(lintegral, cintegral);
2105: PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
2106: return(0);
2107: }
2109: /*@
2110: DMPlexComputeCellwiseIntegralFEM - Form the vector of cellwise integrals F from the global input X using pointwise functions specified by the user
2112: Input Parameters:
2113: + dm - The mesh
2114: . X - Global input vector
2115: - user - The user context
2117: Output Parameter:
2118: . integral - Cellwise integrals for each field
2120: Level: developer
2122: .seealso: DMPlexComputeResidualFEM()
2123: @*/
2124: PetscErrorCode DMPlexComputeCellwiseIntegralFEM(DM dm, Vec X, Vec F, void *user)
2125: {
2126: DM_Plex *mesh = (DM_Plex *) dm->data;
2127: DM dmF;
2128: PetscSection sectionF;
2129: PetscScalar *cintegral, *af;
2130: PetscInt Nf, f, cellHeight, cStart, cEnd, cell;
2137: PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
2138: DMGetNumFields(dm, &Nf);
2139: DMPlexGetVTKCellHeight(dm, &cellHeight);
2140: DMPlexGetSimplexOrBoxCells(dm, cellHeight, &cStart, &cEnd);
2141: /* TODO Introduce a loop over large chunks (right now this is a single chunk) */
2142: PetscCalloc1((cEnd-cStart)*Nf, &cintegral);
2143: DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);
2144: /* Put values in F*/
2145: VecGetDM(F, &dmF);
2146: DMGetLocalSection(dmF, §ionF);
2147: VecGetArray(F, &af);
2148: for (cell = cStart; cell < cEnd; ++cell) {
2149: const PetscInt c = cell - cStart;
2150: PetscInt dof, off;
2152: if (mesh->printFEM > 1) {DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);}
2153: PetscSectionGetDof(sectionF, cell, &dof);
2154: PetscSectionGetOffset(sectionF, cell, &off);
2155: if (dof != Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cell dofs %D != %D", dof, Nf);
2156: for (f = 0; f < Nf; ++f) af[off+f] = cintegral[c*Nf+f];
2157: }
2158: VecRestoreArray(F, &af);
2159: PetscFree(cintegral);
2160: PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
2161: return(0);
2162: }
2164: static PetscErrorCode DMPlexComputeBdIntegral_Internal(DM dm, Vec locX, IS pointIS,
2165: void (*func)(PetscInt, PetscInt, PetscInt,
2166: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
2167: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
2168: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
2169: PetscScalar *fintegral, void *user)
2170: {
2171: DM plex = NULL, plexA = NULL;
2172: DMEnclosureType encAux;
2173: PetscDS prob, probAux = NULL;
2174: PetscSection section, sectionAux = NULL;
2175: Vec locA = NULL;
2176: DMField coordField;
2177: PetscInt Nf, totDim, *uOff, *uOff_x;
2178: PetscInt NfAux = 0, totDimAux = 0, *aOff = NULL;
2179: PetscScalar *u, *a = NULL;
2180: const PetscScalar *constants;
2181: PetscInt numConstants, f;
2182: PetscErrorCode ierr;
2185: DMGetCoordinateField(dm, &coordField);
2186: DMConvert(dm, DMPLEX, &plex);
2187: DMGetDS(dm, &prob);
2188: DMGetLocalSection(dm, §ion);
2189: PetscSectionGetNumFields(section, &Nf);
2190: /* Determine which discretizations we have */
2191: for (f = 0; f < Nf; ++f) {
2192: PetscObject obj;
2193: PetscClassId id;
2195: PetscDSGetDiscretization(prob, f, &obj);
2196: PetscObjectGetClassId(obj, &id);
2197: if (id == PETSCFV_CLASSID) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Not supported for FVM (field %D)", f);
2198: }
2199: /* Read DS information */
2200: PetscDSGetTotalDimension(prob, &totDim);
2201: PetscDSGetComponentOffsets(prob, &uOff);
2202: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
2203: PetscDSGetConstants(prob, &numConstants, &constants);
2204: /* Read Auxiliary DS information */
2205: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
2206: if (locA) {
2207: DM dmAux;
2209: VecGetDM(locA, &dmAux);
2210: DMGetEnclosureRelation(dmAux, dm, &encAux);
2211: DMConvert(dmAux, DMPLEX, &plexA);
2212: DMGetDS(dmAux, &probAux);
2213: PetscDSGetNumFields(probAux, &NfAux);
2214: DMGetLocalSection(dmAux, §ionAux);
2215: PetscDSGetTotalDimension(probAux, &totDimAux);
2216: PetscDSGetComponentOffsets(probAux, &aOff);
2217: }
2218: /* Integrate over points */
2219: {
2220: PetscFEGeom *fgeom, *chunkGeom = NULL;
2221: PetscInt maxDegree;
2222: PetscQuadrature qGeom = NULL;
2223: const PetscInt *points;
2224: PetscInt numFaces, face, Nq, field;
2225: PetscInt numChunks, chunkSize, chunk, Nr, offset;
2227: ISGetLocalSize(pointIS, &numFaces);
2228: ISGetIndices(pointIS, &points);
2229: PetscCalloc2(numFaces*totDim, &u, locA ? numFaces*totDimAux : 0, &a);
2230: DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree);
2231: for (field = 0; field < Nf; ++field) {
2232: PetscFE fe;
2234: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
2235: if (maxDegree <= 1) {DMFieldCreateDefaultQuadrature(coordField, pointIS, &qGeom);}
2236: if (!qGeom) {
2237: PetscFEGetFaceQuadrature(fe, &qGeom);
2238: PetscObjectReference((PetscObject) qGeom);
2239: }
2240: PetscQuadratureGetData(qGeom, NULL, NULL, &Nq, NULL, NULL);
2241: DMPlexGetFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom);
2242: for (face = 0; face < numFaces; ++face) {
2243: const PetscInt point = points[face], *support;
2244: PetscScalar *x = NULL;
2245: PetscInt i;
2247: DMPlexGetSupport(dm, point, &support);
2248: DMPlexVecGetClosure(plex, section, locX, support[0], NULL, &x);
2249: for (i = 0; i < totDim; ++i) u[face*totDim+i] = x[i];
2250: DMPlexVecRestoreClosure(plex, section, locX, support[0], NULL, &x);
2251: if (locA) {
2252: PetscInt subp;
2253: DMGetEnclosurePoint(plexA, dm, encAux, support[0], &subp);
2254: DMPlexVecGetClosure(plexA, sectionAux, locA, subp, NULL, &x);
2255: for (i = 0; i < totDimAux; ++i) a[f*totDimAux+i] = x[i];
2256: DMPlexVecRestoreClosure(plexA, sectionAux, locA, subp, NULL, &x);
2257: }
2258: }
2259: /* Get blocking */
2260: {
2261: PetscQuadrature q;
2262: PetscInt numBatches, batchSize, numBlocks, blockSize;
2263: PetscInt Nq, Nb;
2265: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
2266: PetscFEGetQuadrature(fe, &q);
2267: PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);
2268: PetscFEGetDimension(fe, &Nb);
2269: blockSize = Nb*Nq;
2270: batchSize = numBlocks * blockSize;
2271: chunkSize = numBatches*batchSize;
2272: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2273: numChunks = numFaces / chunkSize;
2274: Nr = numFaces % chunkSize;
2275: offset = numFaces - Nr;
2276: }
2277: /* Do integration for each field */
2278: for (chunk = 0; chunk < numChunks; ++chunk) {
2279: PetscFEGeomGetChunk(fgeom, chunk*chunkSize, (chunk+1)*chunkSize, &chunkGeom);
2280: PetscFEIntegrateBd(prob, field, func, chunkSize, chunkGeom, u, probAux, a, fintegral);
2281: PetscFEGeomRestoreChunk(fgeom, 0, offset, &chunkGeom);
2282: }
2283: PetscFEGeomGetChunk(fgeom, offset, numFaces, &chunkGeom);
2284: PetscFEIntegrateBd(prob, field, func, Nr, chunkGeom, &u[offset*totDim], probAux, a ? &a[offset*totDimAux] : NULL, &fintegral[offset*Nf]);
2285: PetscFEGeomRestoreChunk(fgeom, offset, numFaces, &chunkGeom);
2286: /* Cleanup data arrays */
2287: DMPlexRestoreFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom);
2288: PetscQuadratureDestroy(&qGeom);
2289: PetscFree2(u, a);
2290: ISRestoreIndices(pointIS, &points);
2291: }
2292: }
2293: if (plex) {DMDestroy(&plex);}
2294: if (plexA) {DMDestroy(&plexA);}
2295: return(0);
2296: }
2298: /*@
2299: DMPlexComputeBdIntegral - Form the integral over the specified boundary from the global input X using pointwise functions specified by the user
2301: Input Parameters:
2302: + dm - The mesh
2303: . X - Global input vector
2304: . label - The boundary DMLabel
2305: . numVals - The number of label values to use, or PETSC_DETERMINE for all values
2306: . vals - The label values to use, or PETSC_NULL for all values
2307: . func = The function to integrate along the boundary
2308: - user - The user context
2310: Output Parameter:
2311: . integral - Integral for each field
2313: Level: developer
2315: .seealso: DMPlexComputeIntegralFEM(), DMPlexComputeBdResidualFEM()
2316: @*/
2317: PetscErrorCode DMPlexComputeBdIntegral(DM dm, Vec X, DMLabel label, PetscInt numVals, const PetscInt vals[],
2318: void (*func)(PetscInt, PetscInt, PetscInt,
2319: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
2320: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
2321: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
2322: PetscScalar *integral, void *user)
2323: {
2324: Vec locX;
2325: PetscSection section;
2326: DMLabel depthLabel;
2327: IS facetIS;
2328: PetscInt dim, Nf, f, v;
2337: PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
2338: DMPlexGetDepthLabel(dm, &depthLabel);
2339: DMGetDimension(dm, &dim);
2340: DMLabelGetStratumIS(depthLabel, dim-1, &facetIS);
2341: DMGetLocalSection(dm, §ion);
2342: PetscSectionGetNumFields(section, &Nf);
2343: /* Get local solution with boundary values */
2344: DMGetLocalVector(dm, &locX);
2345: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);
2346: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
2347: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
2348: /* Loop over label values */
2349: PetscArrayzero(integral, Nf);
2350: for (v = 0; v < numVals; ++v) {
2351: IS pointIS;
2352: PetscInt numFaces, face;
2353: PetscScalar *fintegral;
2355: DMLabelGetStratumIS(label, vals[v], &pointIS);
2356: if (!pointIS) continue; /* No points with that id on this process */
2357: {
2358: IS isectIS;
2360: /* TODO: Special cases of ISIntersect where it is quick to check a priori if one is a superset of the other */
2361: ISIntersect_Caching_Internal(facetIS, pointIS, &isectIS);
2362: ISDestroy(&pointIS);
2363: pointIS = isectIS;
2364: }
2365: ISGetLocalSize(pointIS, &numFaces);
2366: PetscCalloc1(numFaces*Nf, &fintegral);
2367: DMPlexComputeBdIntegral_Internal(dm, locX, pointIS, func, fintegral, user);
2368: /* Sum point contributions into integral */
2369: for (f = 0; f < Nf; ++f) for (face = 0; face < numFaces; ++face) integral[f] += fintegral[face*Nf+f];
2370: PetscFree(fintegral);
2371: ISDestroy(&pointIS);
2372: }
2373: DMRestoreLocalVector(dm, &locX);
2374: ISDestroy(&facetIS);
2375: PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
2376: return(0);
2377: }
2379: /*@
2380: DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to a uniformly refined DM.
2382: Input Parameters:
2383: + dmc - The coarse mesh
2384: . dmf - The fine mesh
2385: . isRefined - Flag indicating regular refinement, rather than the same topology
2386: - user - The user context
2388: Output Parameter:
2389: . In - The interpolation matrix
2391: Level: developer
2393: .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
2394: @*/
2395: PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, PetscBool isRefined, Mat In, void *user)
2396: {
2397: DM_Plex *mesh = (DM_Plex *) dmc->data;
2398: const char *name = "Interpolator";
2399: PetscDS cds, rds;
2400: PetscFE *feRef;
2401: PetscFV *fvRef;
2402: PetscSection fsection, fglobalSection;
2403: PetscSection csection, cglobalSection;
2404: PetscScalar *elemMat;
2405: PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, c;
2406: PetscInt cTotDim, rTotDim = 0;
2407: PetscErrorCode ierr;
2410: PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
2411: DMGetDimension(dmf, &dim);
2412: DMGetLocalSection(dmf, &fsection);
2413: DMGetGlobalSection(dmf, &fglobalSection);
2414: DMGetLocalSection(dmc, &csection);
2415: DMGetGlobalSection(dmc, &cglobalSection);
2416: PetscSectionGetNumFields(fsection, &Nf);
2417: DMPlexGetSimplexOrBoxCells(dmc, 0, &cStart, &cEnd);
2418: DMGetDS(dmc, &cds);
2419: DMGetDS(dmf, &rds);
2420: PetscCalloc2(Nf, &feRef, Nf, &fvRef);
2421: for (f = 0; f < Nf; ++f) {
2422: PetscObject obj;
2423: PetscClassId id;
2424: PetscInt rNb = 0, Nc = 0;
2426: PetscDSGetDiscretization(rds, f, &obj);
2427: PetscObjectGetClassId(obj, &id);
2428: if (id == PETSCFE_CLASSID) {
2429: PetscFE fe = (PetscFE) obj;
2431: if (isRefined) {
2432: PetscFERefine(fe, &feRef[f]);
2433: } else {
2434: PetscObjectReference((PetscObject) fe);
2435: feRef[f] = fe;
2436: }
2437: PetscFEGetDimension(feRef[f], &rNb);
2438: PetscFEGetNumComponents(fe, &Nc);
2439: } else if (id == PETSCFV_CLASSID) {
2440: PetscFV fv = (PetscFV) obj;
2441: PetscDualSpace Q;
2443: if (isRefined) {
2444: PetscFVRefine(fv, &fvRef[f]);
2445: } else {
2446: PetscObjectReference((PetscObject) fv);
2447: fvRef[f] = fv;
2448: }
2449: PetscFVGetDualSpace(fvRef[f], &Q);
2450: PetscDualSpaceGetDimension(Q, &rNb);
2451: PetscFVGetNumComponents(fv, &Nc);
2452: }
2453: rTotDim += rNb;
2454: }
2455: PetscDSGetTotalDimension(cds, &cTotDim);
2456: PetscMalloc1(rTotDim*cTotDim,&elemMat);
2457: PetscArrayzero(elemMat, rTotDim*cTotDim);
2458: for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
2459: PetscDualSpace Qref;
2460: PetscQuadrature f;
2461: const PetscReal *qpoints, *qweights;
2462: PetscReal *points;
2463: PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d;
2465: /* Compose points from all dual basis functionals */
2466: if (feRef[fieldI]) {
2467: PetscFEGetDualSpace(feRef[fieldI], &Qref);
2468: PetscFEGetNumComponents(feRef[fieldI], &Nc);
2469: } else {
2470: PetscFVGetDualSpace(fvRef[fieldI], &Qref);
2471: PetscFVGetNumComponents(fvRef[fieldI], &Nc);
2472: }
2473: PetscDualSpaceGetDimension(Qref, &fpdim);
2474: for (i = 0; i < fpdim; ++i) {
2475: PetscDualSpaceGetFunctional(Qref, i, &f);
2476: PetscQuadratureGetData(f, NULL, NULL, &Np, NULL, NULL);
2477: npoints += Np;
2478: }
2479: PetscMalloc1(npoints*dim,&points);
2480: for (i = 0, k = 0; i < fpdim; ++i) {
2481: PetscDualSpaceGetFunctional(Qref, i, &f);
2482: PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
2483: for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
2484: }
2486: for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
2487: PetscObject obj;
2488: PetscClassId id;
2489: PetscInt NcJ = 0, cpdim = 0, j, qNc;
2491: PetscDSGetDiscretization(cds, fieldJ, &obj);
2492: PetscObjectGetClassId(obj, &id);
2493: if (id == PETSCFE_CLASSID) {
2494: PetscFE fe = (PetscFE) obj;
2495: PetscTabulation T = NULL;
2497: /* Evaluate basis at points */
2498: PetscFEGetNumComponents(fe, &NcJ);
2499: PetscFEGetDimension(fe, &cpdim);
2500: /* For now, fields only interpolate themselves */
2501: if (fieldI == fieldJ) {
2502: 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);
2503: PetscFECreateTabulation(fe, 1, npoints, points, 0, &T);
2504: for (i = 0, k = 0; i < fpdim; ++i) {
2505: PetscDualSpaceGetFunctional(Qref, i, &f);
2506: PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
2507: 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);
2508: for (p = 0; p < Np; ++p, ++k) {
2509: for (j = 0; j < cpdim; ++j) {
2510: /*
2511: cTotDim: Total columns in element interpolation matrix, sum of number of dual basis functionals in each field
2512: offsetI, offsetJ: Offsets into the larger element interpolation matrix for different fields
2513: fpdim, i, cpdim, j: Dofs for fine and coarse grids, correspond to dual space basis functionals
2514: qNC, Nc, Ncj, c: Number of components in this field
2515: Np, p: Number of quad points in the fine grid functional i
2516: k: i*Np + p, overall point number for the interpolation
2517: */
2518: for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += T->T[0][k*cpdim*NcJ+j*Nc+c]*qweights[p*qNc+c];
2519: }
2520: }
2521: }
2522: PetscTabulationDestroy(&T);
2523: }
2524: } else if (id == PETSCFV_CLASSID) {
2525: PetscFV fv = (PetscFV) obj;
2527: /* Evaluate constant function at points */
2528: PetscFVGetNumComponents(fv, &NcJ);
2529: cpdim = 1;
2530: /* For now, fields only interpolate themselves */
2531: if (fieldI == fieldJ) {
2532: 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);
2533: for (i = 0, k = 0; i < fpdim; ++i) {
2534: PetscDualSpaceGetFunctional(Qref, i, &f);
2535: PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
2536: 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);
2537: for (p = 0; p < Np; ++p, ++k) {
2538: for (j = 0; j < cpdim; ++j) {
2539: for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += 1.0*qweights[p*qNc+c];
2540: }
2541: }
2542: }
2543: }
2544: }
2545: offsetJ += cpdim;
2546: }
2547: offsetI += fpdim;
2548: PetscFree(points);
2549: }
2550: if (mesh->printFEM > 1) {DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);}
2551: /* Preallocate matrix */
2552: {
2553: Mat preallocator;
2554: PetscScalar *vals;
2555: PetscInt *cellCIndices, *cellFIndices;
2556: PetscInt locRows, locCols, cell;
2558: MatGetLocalSize(In, &locRows, &locCols);
2559: MatCreate(PetscObjectComm((PetscObject) In), &preallocator);
2560: MatSetType(preallocator, MATPREALLOCATOR);
2561: MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);
2562: MatSetUp(preallocator);
2563: PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
2564: for (cell = cStart; cell < cEnd; ++cell) {
2565: if (isRefined) {
2566: DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
2567: MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);
2568: } else {
2569: DMPlexMatSetClosureGeneral(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, preallocator, cell, vals, INSERT_VALUES);
2570: }
2571: }
2572: PetscFree3(vals,cellCIndices,cellFIndices);
2573: MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
2574: MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
2575: MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);
2576: MatDestroy(&preallocator);
2577: }
2578: /* Fill matrix */
2579: MatZeroEntries(In);
2580: for (c = cStart; c < cEnd; ++c) {
2581: if (isRefined) {
2582: DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
2583: } else {
2584: DMPlexMatSetClosureGeneral(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
2585: }
2586: }
2587: for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
2588: PetscFree2(feRef,fvRef);
2589: PetscFree(elemMat);
2590: MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
2591: MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
2592: if (mesh->printFEM) {
2593: PetscPrintf(PetscObjectComm((PetscObject)In), "%s:\n", name);
2594: MatChop(In, 1.0e-10);
2595: MatView(In, NULL);
2596: }
2597: PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
2598: return(0);
2599: }
2601: PetscErrorCode DMPlexComputeMassMatrixNested(DM dmc, DM dmf, Mat mass, void *user)
2602: {
2603: SETERRQ(PetscObjectComm((PetscObject) dmc), PETSC_ERR_SUP, "Laziness");
2604: }
2606: /*@
2607: DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM.
2609: Input Parameters:
2610: + dmf - The fine mesh
2611: . dmc - The coarse mesh
2612: - user - The user context
2614: Output Parameter:
2615: . In - The interpolation matrix
2617: Level: developer
2619: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
2620: @*/
2621: PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user)
2622: {
2623: DM_Plex *mesh = (DM_Plex *) dmf->data;
2624: const char *name = "Interpolator";
2625: PetscDS prob;
2626: PetscSection fsection, csection, globalFSection, globalCSection;
2627: PetscHSetIJ ht;
2628: PetscLayout rLayout;
2629: PetscInt *dnz, *onz;
2630: PetscInt locRows, rStart, rEnd;
2631: PetscReal *x, *v0, *J, *invJ, detJ;
2632: PetscReal *v0c, *Jc, *invJc, detJc;
2633: PetscScalar *elemMat;
2634: PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell;
2638: PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
2639: DMGetCoordinateDim(dmc, &dim);
2640: DMGetDS(dmc, &prob);
2641: PetscDSGetWorkspace(prob, &x, NULL, NULL, NULL, NULL);
2642: PetscDSGetNumFields(prob, &Nf);
2643: PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
2644: PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
2645: DMGetLocalSection(dmf, &fsection);
2646: DMGetGlobalSection(dmf, &globalFSection);
2647: DMGetLocalSection(dmc, &csection);
2648: DMGetGlobalSection(dmc, &globalCSection);
2649: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
2650: PetscDSGetTotalDimension(prob, &totDim);
2651: PetscMalloc1(totDim, &elemMat);
2653: MatGetLocalSize(In, &locRows, NULL);
2654: PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
2655: PetscLayoutSetLocalSize(rLayout, locRows);
2656: PetscLayoutSetBlockSize(rLayout, 1);
2657: PetscLayoutSetUp(rLayout);
2658: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
2659: PetscLayoutDestroy(&rLayout);
2660: PetscCalloc2(locRows,&dnz,locRows,&onz);
2661: PetscHSetIJCreate(&ht);
2662: for (field = 0; field < Nf; ++field) {
2663: PetscObject obj;
2664: PetscClassId id;
2665: PetscDualSpace Q = NULL;
2666: PetscQuadrature f;
2667: const PetscReal *qpoints;
2668: PetscInt Nc, Np, fpdim, i, d;
2670: PetscDSGetDiscretization(prob, field, &obj);
2671: PetscObjectGetClassId(obj, &id);
2672: if (id == PETSCFE_CLASSID) {
2673: PetscFE fe = (PetscFE) obj;
2675: PetscFEGetDualSpace(fe, &Q);
2676: PetscFEGetNumComponents(fe, &Nc);
2677: } else if (id == PETSCFV_CLASSID) {
2678: PetscFV fv = (PetscFV) obj;
2680: PetscFVGetDualSpace(fv, &Q);
2681: Nc = 1;
2682: }
2683: PetscDualSpaceGetDimension(Q, &fpdim);
2684: /* For each fine grid cell */
2685: for (cell = cStart; cell < cEnd; ++cell) {
2686: PetscInt *findices, *cindices;
2687: PetscInt numFIndices, numCIndices;
2689: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
2690: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2691: if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %D != %D dual basis vecs", numFIndices, fpdim);
2692: for (i = 0; i < fpdim; ++i) {
2693: Vec pointVec;
2694: PetscScalar *pV;
2695: PetscSF coarseCellSF = NULL;
2696: const PetscSFNode *coarseCells;
2697: PetscInt numCoarseCells, q, c;
2699: /* Get points from the dual basis functional quadrature */
2700: PetscDualSpaceGetFunctional(Q, i, &f);
2701: PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
2702: VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
2703: VecSetBlockSize(pointVec, dim);
2704: VecGetArray(pointVec, &pV);
2705: for (q = 0; q < Np; ++q) {
2706: const PetscReal xi0[3] = {-1., -1., -1.};
2708: /* Transform point to real space */
2709: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2710: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2711: }
2712: VecRestoreArray(pointVec, &pV);
2713: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2714: /* OPT: Pack all quad points from fine cell */
2715: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2716: PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
2717: /* Update preallocation info */
2718: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2719: if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2720: {
2721: PetscHashIJKey key;
2722: PetscBool missing;
2724: key.i = findices[i];
2725: if (key.i >= 0) {
2726: /* Get indices for coarse elements */
2727: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2728: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);
2729: for (c = 0; c < numCIndices; ++c) {
2730: key.j = cindices[c];
2731: if (key.j < 0) continue;
2732: PetscHSetIJQueryAdd(ht, key, &missing);
2733: if (missing) {
2734: if ((key.j >= rStart) && (key.j < rEnd)) ++dnz[key.i-rStart];
2735: else ++onz[key.i-rStart];
2736: }
2737: }
2738: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);
2739: }
2740: }
2741: }
2742: PetscSFDestroy(&coarseCellSF);
2743: VecDestroy(&pointVec);
2744: }
2745: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
2746: }
2747: }
2748: PetscHSetIJDestroy(&ht);
2749: MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);
2750: MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
2751: PetscFree2(dnz,onz);
2752: for (field = 0; field < Nf; ++field) {
2753: PetscObject obj;
2754: PetscClassId id;
2755: PetscDualSpace Q = NULL;
2756: PetscTabulation T = NULL;
2757: PetscQuadrature f;
2758: const PetscReal *qpoints, *qweights;
2759: PetscInt Nc, qNc, Np, fpdim, i, d;
2761: PetscDSGetDiscretization(prob, field, &obj);
2762: PetscObjectGetClassId(obj, &id);
2763: if (id == PETSCFE_CLASSID) {
2764: PetscFE fe = (PetscFE) obj;
2766: PetscFEGetDualSpace(fe, &Q);
2767: PetscFEGetNumComponents(fe, &Nc);
2768: PetscFECreateTabulation(fe, 1, 1, x, 0, &T);
2769: } else if (id == PETSCFV_CLASSID) {
2770: PetscFV fv = (PetscFV) obj;
2772: PetscFVGetDualSpace(fv, &Q);
2773: Nc = 1;
2774: } else SETERRQ1(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONG,"Unknown discretization type for field %D",field);
2775: PetscDualSpaceGetDimension(Q, &fpdim);
2776: /* For each fine grid cell */
2777: for (cell = cStart; cell < cEnd; ++cell) {
2778: PetscInt *findices, *cindices;
2779: PetscInt numFIndices, numCIndices;
2781: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
2782: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2783: if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %D != %D dual basis vecs", numFIndices, fpdim);
2784: for (i = 0; i < fpdim; ++i) {
2785: Vec pointVec;
2786: PetscScalar *pV;
2787: PetscSF coarseCellSF = NULL;
2788: const PetscSFNode *coarseCells;
2789: PetscInt numCoarseCells, cpdim, q, c, j;
2791: /* Get points from the dual basis functional quadrature */
2792: PetscDualSpaceGetFunctional(Q, i, &f);
2793: PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);
2794: 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);
2795: VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
2796: VecSetBlockSize(pointVec, dim);
2797: VecGetArray(pointVec, &pV);
2798: for (q = 0; q < Np; ++q) {
2799: const PetscReal xi0[3] = {-1., -1., -1.};
2801: /* Transform point to real space */
2802: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2803: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2804: }
2805: VecRestoreArray(pointVec, &pV);
2806: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2807: /* OPT: Read this out from preallocation information */
2808: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2809: /* Update preallocation info */
2810: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2811: if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2812: VecGetArray(pointVec, &pV);
2813: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2814: PetscReal pVReal[3];
2815: const PetscReal xi0[3] = {-1., -1., -1.};
2817: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);
2818: /* Transform points from real space to coarse reference space */
2819: DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);
2820: for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]);
2821: CoordinatesRealToRef(dim, dim, xi0, v0c, invJc, pVReal, x);
2823: if (id == PETSCFE_CLASSID) {
2824: PetscFE fe = (PetscFE) obj;
2826: /* Evaluate coarse basis on contained point */
2827: PetscFEGetDimension(fe, &cpdim);
2828: PetscFEComputeTabulation(fe, 1, x, 0, T);
2829: PetscArrayzero(elemMat, cpdim);
2830: /* Get elemMat entries by multiplying by weight */
2831: for (j = 0; j < cpdim; ++j) {
2832: for (c = 0; c < Nc; ++c) elemMat[j] += T->T[0][j*Nc + c]*qweights[ccell*qNc + c];
2833: }
2834: } else {
2835: cpdim = 1;
2836: for (j = 0; j < cpdim; ++j) {
2837: for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c];
2838: }
2839: }
2840: /* Update interpolator */
2841: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2842: if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2843: MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);
2844: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);
2845: }
2846: VecRestoreArray(pointVec, &pV);
2847: PetscSFDestroy(&coarseCellSF);
2848: VecDestroy(&pointVec);
2849: }
2850: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
2851: }
2852: if (id == PETSCFE_CLASSID) {PetscTabulationDestroy(&T);}
2853: }
2854: PetscFree3(v0,J,invJ);
2855: PetscFree3(v0c,Jc,invJc);
2856: PetscFree(elemMat);
2857: MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
2858: MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
2859: PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
2860: return(0);
2861: }
2863: /*@
2864: DMPlexComputeMassMatrixGeneral - Form the local portion of the mass matrix M from the coarse DM to a non-nested fine DM.
2866: Input Parameters:
2867: + dmf - The fine mesh
2868: . dmc - The coarse mesh
2869: - user - The user context
2871: Output Parameter:
2872: . mass - The mass matrix
2874: Level: developer
2876: .seealso: DMPlexComputeMassMatrixNested(), DMPlexComputeInterpolatorNested(), DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
2877: @*/
2878: PetscErrorCode DMPlexComputeMassMatrixGeneral(DM dmc, DM dmf, Mat mass, void *user)
2879: {
2880: DM_Plex *mesh = (DM_Plex *) dmf->data;
2881: const char *name = "Mass Matrix";
2882: PetscDS prob;
2883: PetscSection fsection, csection, globalFSection, globalCSection;
2884: PetscHSetIJ ht;
2885: PetscLayout rLayout;
2886: PetscInt *dnz, *onz;
2887: PetscInt locRows, rStart, rEnd;
2888: PetscReal *x, *v0, *J, *invJ, detJ;
2889: PetscReal *v0c, *Jc, *invJc, detJc;
2890: PetscScalar *elemMat;
2891: PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell;
2895: DMGetCoordinateDim(dmc, &dim);
2896: DMGetDS(dmc, &prob);
2897: PetscDSGetWorkspace(prob, &x, NULL, NULL, NULL, NULL);
2898: PetscDSGetNumFields(prob, &Nf);
2899: PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
2900: PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
2901: DMGetLocalSection(dmf, &fsection);
2902: DMGetGlobalSection(dmf, &globalFSection);
2903: DMGetLocalSection(dmc, &csection);
2904: DMGetGlobalSection(dmc, &globalCSection);
2905: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
2906: PetscDSGetTotalDimension(prob, &totDim);
2907: PetscMalloc1(totDim, &elemMat);
2909: MatGetLocalSize(mass, &locRows, NULL);
2910: PetscLayoutCreate(PetscObjectComm((PetscObject) mass), &rLayout);
2911: PetscLayoutSetLocalSize(rLayout, locRows);
2912: PetscLayoutSetBlockSize(rLayout, 1);
2913: PetscLayoutSetUp(rLayout);
2914: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
2915: PetscLayoutDestroy(&rLayout);
2916: PetscCalloc2(locRows,&dnz,locRows,&onz);
2917: PetscHSetIJCreate(&ht);
2918: for (field = 0; field < Nf; ++field) {
2919: PetscObject obj;
2920: PetscClassId id;
2921: PetscQuadrature quad;
2922: const PetscReal *qpoints;
2923: PetscInt Nq, Nc, i, d;
2925: PetscDSGetDiscretization(prob, field, &obj);
2926: PetscObjectGetClassId(obj, &id);
2927: if (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, &quad);}
2928: else {PetscFVGetQuadrature((PetscFV) obj, &quad);}
2929: PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, NULL);
2930: /* For each fine grid cell */
2931: for (cell = cStart; cell < cEnd; ++cell) {
2932: Vec pointVec;
2933: PetscScalar *pV;
2934: PetscSF coarseCellSF = NULL;
2935: const PetscSFNode *coarseCells;
2936: PetscInt numCoarseCells, q, c;
2937: PetscInt *findices, *cindices;
2938: PetscInt numFIndices, numCIndices;
2940: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
2941: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2942: /* Get points from the quadrature */
2943: VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);
2944: VecSetBlockSize(pointVec, dim);
2945: VecGetArray(pointVec, &pV);
2946: for (q = 0; q < Nq; ++q) {
2947: const PetscReal xi0[3] = {-1., -1., -1.};
2949: /* Transform point to real space */
2950: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2951: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2952: }
2953: VecRestoreArray(pointVec, &pV);
2954: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2955: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2956: PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
2957: /* Update preallocation info */
2958: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2959: if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2960: {
2961: PetscHashIJKey key;
2962: PetscBool missing;
2964: for (i = 0; i < numFIndices; ++i) {
2965: key.i = findices[i];
2966: if (key.i >= 0) {
2967: /* Get indices for coarse elements */
2968: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2969: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);
2970: for (c = 0; c < numCIndices; ++c) {
2971: key.j = cindices[c];
2972: if (key.j < 0) continue;
2973: PetscHSetIJQueryAdd(ht, key, &missing);
2974: if (missing) {
2975: if ((key.j >= rStart) && (key.j < rEnd)) ++dnz[key.i-rStart];
2976: else ++onz[key.i-rStart];
2977: }
2978: }
2979: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);
2980: }
2981: }
2982: }
2983: }
2984: PetscSFDestroy(&coarseCellSF);
2985: VecDestroy(&pointVec);
2986: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
2987: }
2988: }
2989: PetscHSetIJDestroy(&ht);
2990: MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
2991: MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
2992: PetscFree2(dnz,onz);
2993: for (field = 0; field < Nf; ++field) {
2994: PetscObject obj;
2995: PetscClassId id;
2996: PetscTabulation T, Tfine;
2997: PetscQuadrature quad;
2998: const PetscReal *qpoints, *qweights;
2999: PetscInt Nq, Nc, i, d;
3001: PetscDSGetDiscretization(prob, field, &obj);
3002: PetscObjectGetClassId(obj, &id);
3003: if (id == PETSCFE_CLASSID) {
3004: PetscFEGetQuadrature((PetscFE) obj, &quad);
3005: PetscFEGetCellTabulation((PetscFE) obj, &Tfine);
3006: PetscFECreateTabulation((PetscFE) obj, 1, 1, x, 0, &T);
3007: } else {
3008: PetscFVGetQuadrature((PetscFV) obj, &quad);
3009: }
3010: PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, &qweights);
3011: /* For each fine grid cell */
3012: for (cell = cStart; cell < cEnd; ++cell) {
3013: Vec pointVec;
3014: PetscScalar *pV;
3015: PetscSF coarseCellSF = NULL;
3016: const PetscSFNode *coarseCells;
3017: PetscInt numCoarseCells, cpdim, q, c, j;
3018: PetscInt *findices, *cindices;
3019: PetscInt numFIndices, numCIndices;
3021: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
3022: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
3023: /* Get points from the quadrature */
3024: VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);
3025: VecSetBlockSize(pointVec, dim);
3026: VecGetArray(pointVec, &pV);
3027: for (q = 0; q < Nq; ++q) {
3028: const PetscReal xi0[3] = {-1., -1., -1.};
3030: /* Transform point to real space */
3031: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
3032: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
3033: }
3034: VecRestoreArray(pointVec, &pV);
3035: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
3036: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
3037: /* Update matrix */
3038: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
3039: if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
3040: VecGetArray(pointVec, &pV);
3041: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
3042: PetscReal pVReal[3];
3043: const PetscReal xi0[3] = {-1., -1., -1.};
3046: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);
3047: /* Transform points from real space to coarse reference space */
3048: DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);
3049: for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]);
3050: CoordinatesRealToRef(dim, dim, xi0, v0c, invJc, pVReal, x);
3052: if (id == PETSCFE_CLASSID) {
3053: PetscFE fe = (PetscFE) obj;
3055: /* Evaluate coarse basis on contained point */
3056: PetscFEGetDimension(fe, &cpdim);
3057: PetscFEComputeTabulation(fe, 1, x, 0, T);
3058: /* Get elemMat entries by multiplying by weight */
3059: for (i = 0; i < numFIndices; ++i) {
3060: PetscArrayzero(elemMat, cpdim);
3061: for (j = 0; j < cpdim; ++j) {
3062: for (c = 0; c < Nc; ++c) elemMat[j] += T->T[0][j*Nc + c]*Tfine->T[0][(ccell*numFIndices + i)*Nc + c]*qweights[ccell*Nc + c]*detJ;
3063: }
3064: /* Update interpolator */
3065: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
3066: if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
3067: MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);
3068: }
3069: } else {
3070: cpdim = 1;
3071: for (i = 0; i < numFIndices; ++i) {
3072: PetscArrayzero(elemMat, cpdim);
3073: for (j = 0; j < cpdim; ++j) {
3074: for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*1.0*qweights[ccell*Nc + c]*detJ;
3075: }
3076: /* Update interpolator */
3077: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
3078: PetscPrintf(PETSC_COMM_SELF, "Nq: %D %D Nf: %D %D Nc: %D %D\n", ccell, Nq, i, numFIndices, j, numCIndices);
3079: if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
3080: MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);
3081: }
3082: }
3083: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);
3084: }
3085: VecRestoreArray(pointVec, &pV);
3086: PetscSFDestroy(&coarseCellSF);
3087: VecDestroy(&pointVec);
3088: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
3089: }
3090: if (id == PETSCFE_CLASSID) {PetscTabulationDestroy(&T);}
3091: }
3092: PetscFree3(v0,J,invJ);
3093: PetscFree3(v0c,Jc,invJc);
3094: PetscFree(elemMat);
3095: MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
3096: MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
3097: return(0);
3098: }
3100: /*@
3101: DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns
3103: Input Parameters:
3104: + dmc - The coarse mesh
3105: - dmf - The fine mesh
3106: - user - The user context
3108: Output Parameter:
3109: . sc - The mapping
3111: Level: developer
3113: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
3114: @*/
3115: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
3116: {
3117: PetscDS prob;
3118: PetscFE *feRef;
3119: PetscFV *fvRef;
3120: Vec fv, cv;
3121: IS fis, cis;
3122: PetscSection fsection, fglobalSection, csection, cglobalSection;
3123: PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
3124: PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, c, dim, d, startC, endC, offsetC, offsetF, m;
3125: PetscBool *needAvg;
3129: PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
3130: DMGetDimension(dmf, &dim);
3131: DMGetLocalSection(dmf, &fsection);
3132: DMGetGlobalSection(dmf, &fglobalSection);
3133: DMGetLocalSection(dmc, &csection);
3134: DMGetGlobalSection(dmc, &cglobalSection);
3135: PetscSectionGetNumFields(fsection, &Nf);
3136: DMPlexGetSimplexOrBoxCells(dmc, 0, &cStart, &cEnd);
3137: DMGetDS(dmc, &prob);
3138: PetscCalloc3(Nf,&feRef,Nf,&fvRef,Nf,&needAvg);
3139: for (f = 0; f < Nf; ++f) {
3140: PetscObject obj;
3141: PetscClassId id;
3142: PetscInt fNb = 0, Nc = 0;
3144: PetscDSGetDiscretization(prob, f, &obj);
3145: PetscObjectGetClassId(obj, &id);
3146: if (id == PETSCFE_CLASSID) {
3147: PetscFE fe = (PetscFE) obj;
3148: PetscSpace sp;
3149: PetscInt maxDegree;
3151: PetscFERefine(fe, &feRef[f]);
3152: PetscFEGetDimension(feRef[f], &fNb);
3153: PetscFEGetNumComponents(fe, &Nc);
3154: PetscFEGetBasisSpace(fe, &sp);
3155: PetscSpaceGetDegree(sp, NULL, &maxDegree);
3156: if (!maxDegree) needAvg[f] = PETSC_TRUE;
3157: } else if (id == PETSCFV_CLASSID) {
3158: PetscFV fv = (PetscFV) obj;
3159: PetscDualSpace Q;
3161: PetscFVRefine(fv, &fvRef[f]);
3162: PetscFVGetDualSpace(fvRef[f], &Q);
3163: PetscDualSpaceGetDimension(Q, &fNb);
3164: PetscFVGetNumComponents(fv, &Nc);
3165: needAvg[f] = PETSC_TRUE;
3166: }
3167: fTotDim += fNb;
3168: }
3169: PetscDSGetTotalDimension(prob, &cTotDim);
3170: PetscMalloc1(cTotDim,&cmap);
3171: for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
3172: PetscFE feC;
3173: PetscFV fvC;
3174: PetscDualSpace QF, QC;
3175: PetscInt order = -1, NcF, NcC, fpdim, cpdim;
3177: if (feRef[field]) {
3178: PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
3179: PetscFEGetNumComponents(feC, &NcC);
3180: PetscFEGetNumComponents(feRef[field], &NcF);
3181: PetscFEGetDualSpace(feRef[field], &QF);
3182: PetscDualSpaceGetOrder(QF, &order);
3183: PetscDualSpaceGetDimension(QF, &fpdim);
3184: PetscFEGetDualSpace(feC, &QC);
3185: PetscDualSpaceGetDimension(QC, &cpdim);
3186: } else {
3187: PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);
3188: PetscFVGetNumComponents(fvC, &NcC);
3189: PetscFVGetNumComponents(fvRef[field], &NcF);
3190: PetscFVGetDualSpace(fvRef[field], &QF);
3191: PetscDualSpaceGetDimension(QF, &fpdim);
3192: PetscFVGetDualSpace(fvC, &QC);
3193: PetscDualSpaceGetDimension(QC, &cpdim);
3194: }
3195: 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);
3196: for (c = 0; c < cpdim; ++c) {
3197: PetscQuadrature cfunc;
3198: const PetscReal *cqpoints, *cqweights;
3199: PetscInt NqcC, NpC;
3200: PetscBool found = PETSC_FALSE;
3202: PetscDualSpaceGetFunctional(QC, c, &cfunc);
3203: PetscQuadratureGetData(cfunc, NULL, &NqcC, &NpC, &cqpoints, &cqweights);
3204: 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);
3205: if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
3206: for (f = 0; f < fpdim; ++f) {
3207: PetscQuadrature ffunc;
3208: const PetscReal *fqpoints, *fqweights;
3209: PetscReal sum = 0.0;
3210: PetscInt NqcF, NpF;
3212: PetscDualSpaceGetFunctional(QF, f, &ffunc);
3213: PetscQuadratureGetData(ffunc, NULL, &NqcF, &NpF, &fqpoints, &fqweights);
3214: 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);
3215: if (NpC != NpF) continue;
3216: for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
3217: if (sum > 1.0e-9) continue;
3218: for (d = 0; d < NcC; ++d) sum += PetscAbsReal(cqweights[d]*fqweights[d]);
3219: if (sum < 1.0e-9) continue;
3220: cmap[offsetC+c] = offsetF+f;
3221: found = PETSC_TRUE;
3222: break;
3223: }
3224: if (!found) {
3225: /* TODO We really want the average here, but some asshole put VecScatter in the interface */
3226: if (fvRef[field] || (feRef[field] && order == 0)) {
3227: cmap[offsetC+c] = offsetF+0;
3228: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection");
3229: }
3230: }
3231: offsetC += cpdim;
3232: offsetF += fpdim;
3233: }
3234: for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);PetscFVDestroy(&fvRef[f]);}
3235: PetscFree3(feRef,fvRef,needAvg);
3237: DMGetGlobalVector(dmf, &fv);
3238: DMGetGlobalVector(dmc, &cv);
3239: VecGetOwnershipRange(cv, &startC, &endC);
3240: PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
3241: PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
3242: PetscMalloc1(m,&cindices);
3243: PetscMalloc1(m,&findices);
3244: for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
3245: for (c = cStart; c < cEnd; ++c) {
3246: DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
3247: for (d = 0; d < cTotDim; ++d) {
3248: if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue;
3249: 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]]);
3250: cindices[cellCIndices[d]-startC] = cellCIndices[d];
3251: findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
3252: }
3253: }
3254: PetscFree(cmap);
3255: PetscFree2(cellCIndices,cellFIndices);
3257: ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
3258: ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
3259: VecScatterCreate(cv, cis, fv, fis, sc);
3260: ISDestroy(&cis);
3261: ISDestroy(&fis);
3262: DMRestoreGlobalVector(dmf, &fv);
3263: DMRestoreGlobalVector(dmc, &cv);
3264: PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
3265: return(0);
3266: }
3268: /*@C
3269: DMPlexGetCellFields - Retrieve the field values values for a chunk of cells
3271: Input Parameters:
3272: + dm - The DM
3273: . cellIS - The cells to include
3274: . locX - A local vector with the solution fields
3275: . locX_t - A local vector with solution field time derivatives, or NULL
3276: - locA - A local vector with auxiliary fields, or NULL
3278: Output Parameters:
3279: + u - The field coefficients
3280: . u_t - The fields derivative coefficients
3281: - a - The auxiliary field coefficients
3283: Level: developer
3285: .seealso: DMPlexGetFaceFields()
3286: @*/
3287: PetscErrorCode DMPlexGetCellFields(DM dm, IS cellIS, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
3288: {
3289: DM plex, plexA = NULL;
3290: DMEnclosureType encAux;
3291: PetscSection section, sectionAux;
3292: PetscDS prob;
3293: const PetscInt *cells;
3294: PetscInt cStart, cEnd, numCells, totDim, totDimAux, c;
3295: PetscErrorCode ierr;
3305: DMPlexConvertPlex(dm, &plex, PETSC_FALSE);
3306: ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
3307: DMGetLocalSection(dm, §ion);
3308: DMGetCellDS(dm, cells ? cells[cStart] : cStart, &prob);
3309: PetscDSGetTotalDimension(prob, &totDim);
3310: if (locA) {
3311: DM dmAux;
3312: PetscDS probAux;
3314: VecGetDM(locA, &dmAux);
3315: DMGetEnclosureRelation(dmAux, dm, &encAux);
3316: DMPlexConvertPlex(dmAux, &plexA, PETSC_FALSE);
3317: DMGetLocalSection(dmAux, §ionAux);
3318: DMGetDS(dmAux, &probAux);
3319: PetscDSGetTotalDimension(probAux, &totDimAux);
3320: }
3321: numCells = cEnd - cStart;
3322: DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, u);
3323: if (locX_t) {DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, u_t);} else {*u_t = NULL;}
3324: if (locA) {DMGetWorkArray(dm, numCells*totDimAux, MPIU_SCALAR, a);} else {*a = NULL;}
3325: for (c = cStart; c < cEnd; ++c) {
3326: const PetscInt cell = cells ? cells[c] : c;
3327: const PetscInt cind = c - cStart;
3328: PetscScalar *x = NULL, *x_t = NULL, *ul = *u, *ul_t = *u_t, *al = *a;
3329: PetscInt i;
3331: DMPlexVecGetClosure(plex, section, locX, cell, NULL, &x);
3332: for (i = 0; i < totDim; ++i) ul[cind*totDim+i] = x[i];
3333: DMPlexVecRestoreClosure(plex, section, locX, cell, NULL, &x);
3334: if (locX_t) {
3335: DMPlexVecGetClosure(plex, section, locX_t, cell, NULL, &x_t);
3336: for (i = 0; i < totDim; ++i) ul_t[cind*totDim+i] = x_t[i];
3337: DMPlexVecRestoreClosure(plex, section, locX_t, cell, NULL, &x_t);
3338: }
3339: if (locA) {
3340: PetscInt subcell;
3341: DMGetEnclosurePoint(plexA, dm, encAux, cell, &subcell);
3342: DMPlexVecGetClosure(plexA, sectionAux, locA, subcell, NULL, &x);
3343: for (i = 0; i < totDimAux; ++i) al[cind*totDimAux+i] = x[i];
3344: DMPlexVecRestoreClosure(plexA, sectionAux, locA, subcell, NULL, &x);
3345: }
3346: }
3347: DMDestroy(&plex);
3348: if (locA) {DMDestroy(&plexA);}
3349: ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);
3350: return(0);
3351: }
3353: /*@C
3354: DMPlexRestoreCellFields - Restore the field values values for a chunk of cells
3356: Input Parameters:
3357: + dm - The DM
3358: . cellIS - The cells to include
3359: . locX - A local vector with the solution fields
3360: . locX_t - A local vector with solution field time derivatives, or NULL
3361: - locA - A local vector with auxiliary fields, or NULL
3363: Output Parameters:
3364: + u - The field coefficients
3365: . u_t - The fields derivative coefficients
3366: - a - The auxiliary field coefficients
3368: Level: developer
3370: .seealso: DMPlexGetFaceFields()
3371: @*/
3372: PetscErrorCode DMPlexRestoreCellFields(DM dm, IS cellIS, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
3373: {
3377: DMRestoreWorkArray(dm, 0, MPIU_SCALAR, u);
3378: if (locX_t) {DMRestoreWorkArray(dm, 0, MPIU_SCALAR, u_t);}
3379: if (locA) {DMRestoreWorkArray(dm, 0, MPIU_SCALAR, a);}
3380: return(0);
3381: }
3383: /*@C
3384: DMPlexGetFaceFields - Retrieve the field values values for a chunk of faces
3386: Input Parameters:
3387: + dm - The DM
3388: . fStart - The first face to include
3389: . fEnd - The first face to exclude
3390: . locX - A local vector with the solution fields
3391: . locX_t - A local vector with solution field time derivatives, or NULL
3392: . faceGeometry - A local vector with face geometry
3393: . cellGeometry - A local vector with cell geometry
3394: - locaGrad - A local vector with field gradients, or NULL
3396: Output Parameters:
3397: + Nface - The number of faces with field values
3398: . uL - The field values at the left side of the face
3399: - uR - The field values at the right side of the face
3401: Level: developer
3403: .seealso: DMPlexGetCellFields()
3404: @*/
3405: 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)
3406: {
3407: DM dmFace, dmCell, dmGrad = NULL;
3408: PetscSection section;
3409: PetscDS prob;
3410: DMLabel ghostLabel;
3411: const PetscScalar *facegeom, *cellgeom, *x, *lgrad;
3412: PetscBool *isFE;
3413: PetscInt dim, Nf, f, Nc, numFaces = fEnd - fStart, iface, face;
3414: PetscErrorCode ierr;
3425: DMGetDimension(dm, &dim);
3426: DMGetDS(dm, &prob);
3427: DMGetLocalSection(dm, §ion);
3428: PetscDSGetNumFields(prob, &Nf);
3429: PetscDSGetTotalComponents(prob, &Nc);
3430: PetscMalloc1(Nf, &isFE);
3431: for (f = 0; f < Nf; ++f) {
3432: PetscObject obj;
3433: PetscClassId id;
3435: PetscDSGetDiscretization(prob, f, &obj);
3436: PetscObjectGetClassId(obj, &id);
3437: if (id == PETSCFE_CLASSID) {isFE[f] = PETSC_TRUE;}
3438: else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
3439: else {isFE[f] = PETSC_FALSE;}
3440: }
3441: DMGetLabel(dm, "ghost", &ghostLabel);
3442: VecGetArrayRead(locX, &x);
3443: VecGetDM(faceGeometry, &dmFace);
3444: VecGetArrayRead(faceGeometry, &facegeom);
3445: VecGetDM(cellGeometry, &dmCell);
3446: VecGetArrayRead(cellGeometry, &cellgeom);
3447: if (locGrad) {
3448: VecGetDM(locGrad, &dmGrad);
3449: VecGetArrayRead(locGrad, &lgrad);
3450: }
3451: DMGetWorkArray(dm, numFaces*Nc, MPIU_SCALAR, uL);
3452: DMGetWorkArray(dm, numFaces*Nc, MPIU_SCALAR, uR);
3453: /* Right now just eat the extra work for FE (could make a cell loop) */
3454: for (face = fStart, iface = 0; face < fEnd; ++face) {
3455: const PetscInt *cells;
3456: PetscFVFaceGeom *fg;
3457: PetscFVCellGeom *cgL, *cgR;
3458: PetscScalar *xL, *xR, *gL, *gR;
3459: PetscScalar *uLl = *uL, *uRl = *uR;
3460: PetscInt ghost, nsupp, nchild;
3462: DMLabelGetValue(ghostLabel, face, &ghost);
3463: DMPlexGetSupportSize(dm, face, &nsupp);
3464: DMPlexGetTreeChildren(dm, face, &nchild, NULL);
3465: if (ghost >= 0 || nsupp > 2 || nchild > 0) continue;
3466: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
3467: DMPlexGetSupport(dm, face, &cells);
3468: DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
3469: DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
3470: for (f = 0; f < Nf; ++f) {
3471: PetscInt off;
3473: PetscDSGetComponentOffset(prob, f, &off);
3474: if (isFE[f]) {
3475: const PetscInt *cone;
3476: PetscInt comp, coneSizeL, coneSizeR, faceLocL, faceLocR, ldof, rdof, d;
3478: xL = xR = NULL;
3479: PetscSectionGetFieldComponents(section, f, &comp);
3480: DMPlexVecGetClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
3481: DMPlexVecGetClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
3482: DMPlexGetCone(dm, cells[0], &cone);
3483: DMPlexGetConeSize(dm, cells[0], &coneSizeL);
3484: for (faceLocL = 0; faceLocL < coneSizeL; ++faceLocL) if (cone[faceLocL] == face) break;
3485: DMPlexGetCone(dm, cells[1], &cone);
3486: DMPlexGetConeSize(dm, cells[1], &coneSizeR);
3487: for (faceLocR = 0; faceLocR < coneSizeR; ++faceLocR) if (cone[faceLocR] == face) break;
3488: 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]);
3489: /* Check that FEM field has values in the right cell (sometimes its an FV ghost cell) */
3490: /* TODO: this is a hack that might not be right for nonconforming */
3491: if (faceLocL < coneSizeL) {
3492: PetscFEEvaluateFaceFields_Internal(prob, f, faceLocL, xL, &uLl[iface*Nc+off]);
3493: if (rdof == ldof && faceLocR < coneSizeR) {PetscFEEvaluateFaceFields_Internal(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);}
3494: else {for (d = 0; d < comp; ++d) uRl[iface*Nc+off+d] = uLl[iface*Nc+off+d];}
3495: }
3496: else {
3497: PetscFEEvaluateFaceFields_Internal(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);
3498: PetscSectionGetFieldComponents(section, f, &comp);
3499: for (d = 0; d < comp; ++d) uLl[iface*Nc+off+d] = uRl[iface*Nc+off+d];
3500: }
3501: DMPlexVecRestoreClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
3502: DMPlexVecRestoreClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
3503: } else {
3504: PetscFV fv;
3505: PetscInt numComp, c;
3507: PetscDSGetDiscretization(prob, f, (PetscObject *) &fv);
3508: PetscFVGetNumComponents(fv, &numComp);
3509: DMPlexPointLocalFieldRead(dm, cells[0], f, x, &xL);
3510: DMPlexPointLocalFieldRead(dm, cells[1], f, x, &xR);
3511: if (dmGrad) {
3512: PetscReal dxL[3], dxR[3];
3514: DMPlexPointLocalRead(dmGrad, cells[0], lgrad, &gL);
3515: DMPlexPointLocalRead(dmGrad, cells[1], lgrad, &gR);
3516: DMPlex_WaxpyD_Internal(dim, -1, cgL->centroid, fg->centroid, dxL);
3517: DMPlex_WaxpyD_Internal(dim, -1, cgR->centroid, fg->centroid, dxR);
3518: for (c = 0; c < numComp; ++c) {
3519: uLl[iface*Nc+off+c] = xL[c] + DMPlex_DotD_Internal(dim, &gL[c*dim], dxL);
3520: uRl[iface*Nc+off+c] = xR[c] + DMPlex_DotD_Internal(dim, &gR[c*dim], dxR);
3521: }
3522: } else {
3523: for (c = 0; c < numComp; ++c) {
3524: uLl[iface*Nc+off+c] = xL[c];
3525: uRl[iface*Nc+off+c] = xR[c];
3526: }
3527: }
3528: }
3529: }
3530: ++iface;
3531: }
3532: *Nface = iface;
3533: VecRestoreArrayRead(locX, &x);
3534: VecRestoreArrayRead(faceGeometry, &facegeom);
3535: VecRestoreArrayRead(cellGeometry, &cellgeom);
3536: if (locGrad) {
3537: VecRestoreArrayRead(locGrad, &lgrad);
3538: }
3539: PetscFree(isFE);
3540: return(0);
3541: }
3543: /*@C
3544: DMPlexRestoreFaceFields - Restore the field values values for a chunk of faces
3546: Input Parameters:
3547: + dm - The DM
3548: . fStart - The first face to include
3549: . fEnd - The first face to exclude
3550: . locX - A local vector with the solution fields
3551: . locX_t - A local vector with solution field time derivatives, or NULL
3552: . faceGeometry - A local vector with face geometry
3553: . cellGeometry - A local vector with cell geometry
3554: - locaGrad - A local vector with field gradients, or NULL
3556: Output Parameters:
3557: + Nface - The number of faces with field values
3558: . uL - The field values at the left side of the face
3559: - uR - The field values at the right side of the face
3561: Level: developer
3563: .seealso: DMPlexGetFaceFields()
3564: @*/
3565: 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)
3566: {
3570: DMRestoreWorkArray(dm, 0, MPIU_SCALAR, uL);
3571: DMRestoreWorkArray(dm, 0, MPIU_SCALAR, uR);
3572: return(0);
3573: }
3575: /*@C
3576: DMPlexGetFaceGeometry - Retrieve the geometric values for a chunk of faces
3578: Input Parameters:
3579: + dm - The DM
3580: . fStart - The first face to include
3581: . fEnd - The first face to exclude
3582: . faceGeometry - A local vector with face geometry
3583: - cellGeometry - A local vector with cell geometry
3585: Output Parameters:
3586: + Nface - The number of faces with field values
3587: . fgeom - The extract the face centroid and normal
3588: - vol - The cell volume
3590: Level: developer
3592: .seealso: DMPlexGetCellFields()
3593: @*/
3594: PetscErrorCode DMPlexGetFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscInt *Nface, PetscFVFaceGeom **fgeom, PetscReal **vol)
3595: {
3596: DM dmFace, dmCell;
3597: DMLabel ghostLabel;
3598: const PetscScalar *facegeom, *cellgeom;
3599: PetscInt dim, numFaces = fEnd - fStart, iface, face;
3600: PetscErrorCode ierr;
3608: DMGetDimension(dm, &dim);
3609: DMGetLabel(dm, "ghost", &ghostLabel);
3610: VecGetDM(faceGeometry, &dmFace);
3611: VecGetArrayRead(faceGeometry, &facegeom);
3612: VecGetDM(cellGeometry, &dmCell);
3613: VecGetArrayRead(cellGeometry, &cellgeom);
3614: PetscMalloc1(numFaces, fgeom);
3615: DMGetWorkArray(dm, numFaces*2, MPIU_SCALAR, vol);
3616: for (face = fStart, iface = 0; face < fEnd; ++face) {
3617: const PetscInt *cells;
3618: PetscFVFaceGeom *fg;
3619: PetscFVCellGeom *cgL, *cgR;
3620: PetscFVFaceGeom *fgeoml = *fgeom;
3621: PetscReal *voll = *vol;
3622: PetscInt ghost, d, nchild, nsupp;
3624: DMLabelGetValue(ghostLabel, face, &ghost);
3625: DMPlexGetSupportSize(dm, face, &nsupp);
3626: DMPlexGetTreeChildren(dm, face, &nchild, NULL);
3627: if (ghost >= 0 || nsupp > 2 || nchild > 0) continue;
3628: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
3629: DMPlexGetSupport(dm, face, &cells);
3630: DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
3631: DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
3632: for (d = 0; d < dim; ++d) {
3633: fgeoml[iface].centroid[d] = fg->centroid[d];
3634: fgeoml[iface].normal[d] = fg->normal[d];
3635: }
3636: voll[iface*2+0] = cgL->volume;
3637: voll[iface*2+1] = cgR->volume;
3638: ++iface;
3639: }
3640: *Nface = iface;
3641: VecRestoreArrayRead(faceGeometry, &facegeom);
3642: VecRestoreArrayRead(cellGeometry, &cellgeom);
3643: return(0);
3644: }
3646: /*@C
3647: DMPlexRestoreFaceGeometry - Restore the field values values for a chunk of faces
3649: Input Parameters:
3650: + dm - The DM
3651: . fStart - The first face to include
3652: . fEnd - The first face to exclude
3653: . faceGeometry - A local vector with face geometry
3654: - cellGeometry - A local vector with cell geometry
3656: Output Parameters:
3657: + Nface - The number of faces with field values
3658: . fgeom - The extract the face centroid and normal
3659: - vol - The cell volume
3661: Level: developer
3663: .seealso: DMPlexGetFaceFields()
3664: @*/
3665: PetscErrorCode DMPlexRestoreFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscInt *Nface, PetscFVFaceGeom **fgeom, PetscReal **vol)
3666: {
3670: PetscFree(*fgeom);
3671: DMRestoreWorkArray(dm, 0, MPIU_REAL, vol);
3672: return(0);
3673: }
3675: PetscErrorCode DMSNESGetFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
3676: {
3677: char composeStr[33] = {0};
3678: PetscObjectId id;
3679: PetscContainer container;
3680: PetscErrorCode ierr;
3683: PetscObjectGetId((PetscObject)quad,&id);
3684: PetscSNPrintf(composeStr, 32, "DMSNESGetFEGeom_%x\n", id);
3685: PetscObjectQuery((PetscObject) pointIS, composeStr, (PetscObject *) &container);
3686: if (container) {
3687: PetscContainerGetPointer(container, (void **) geom);
3688: } else {
3689: DMFieldCreateFEGeom(coordField, pointIS, quad, faceData, geom);
3690: PetscContainerCreate(PETSC_COMM_SELF,&container);
3691: PetscContainerSetPointer(container, (void *) *geom);
3692: PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom);
3693: PetscObjectCompose((PetscObject) pointIS, composeStr, (PetscObject) container);
3694: PetscContainerDestroy(&container);
3695: }
3696: return(0);
3697: }
3699: PetscErrorCode DMSNESRestoreFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
3700: {
3702: *geom = NULL;
3703: return(0);
3704: }
3706: PetscErrorCode DMPlexComputeResidual_Patch_Internal(DM dm, PetscSection section, IS cellIS, PetscReal t, Vec locX, Vec locX_t, Vec locF, void *user)
3707: {
3708: DM_Plex *mesh = (DM_Plex *) dm->data;
3709: const char *name = "Residual";
3710: DM dmAux = NULL;
3711: DMLabel ghostLabel = NULL;
3712: PetscDS prob = NULL;
3713: PetscDS probAux = NULL;
3714: PetscBool useFEM = PETSC_FALSE;
3715: PetscBool isImplicit = (locX_t || t == PETSC_MIN_REAL) ? PETSC_TRUE : PETSC_FALSE;
3716: DMField coordField = NULL;
3717: Vec locA;
3718: PetscScalar *u = NULL, *u_t, *a, *uL = NULL, *uR = NULL;
3719: IS chunkIS;
3720: const PetscInt *cells;
3721: PetscInt cStart, cEnd, numCells;
3722: PetscInt Nf, f, totDim, totDimAux, numChunks, cellChunkSize, chunk, fStart, fEnd;
3723: PetscInt maxDegree = PETSC_MAX_INT;
3724: PetscQuadrature affineQuad = NULL, *quads = NULL;
3725: PetscFEGeom *affineGeom = NULL, **geoms = NULL;
3726: PetscErrorCode ierr;
3729: PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);
3730: /* FEM+FVM */
3731: /* 1: Get sizes from dm and dmAux */
3732: DMGetLabel(dm, "ghost", &ghostLabel);
3733: DMGetDS(dm, &prob);
3734: PetscDSGetNumFields(prob, &Nf);
3735: PetscDSGetTotalDimension(prob, &totDim);
3736: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
3737: if (locA) {
3738: VecGetDM(locA, &dmAux);
3739: DMGetDS(dmAux, &probAux);
3740: PetscDSGetTotalDimension(probAux, &totDimAux);
3741: }
3742: /* 2: Get geometric data */
3743: for (f = 0; f < Nf; ++f) {
3744: PetscObject obj;
3745: PetscClassId id;
3746: PetscBool fimp;
3748: PetscDSGetImplicit(prob, f, &fimp);
3749: if (isImplicit != fimp) continue;
3750: PetscDSGetDiscretization(prob, f, &obj);
3751: PetscObjectGetClassId(obj, &id);
3752: if (id == PETSCFE_CLASSID) {useFEM = PETSC_TRUE;}
3753: if (id == PETSCFV_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Use of FVM with PCPATCH not yet implemented");
3754: }
3755: if (useFEM) {
3756: DMGetCoordinateField(dm, &coordField);
3757: DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);
3758: if (maxDegree <= 1) {
3759: DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);
3760: if (affineQuad) {
3761: DMSNESGetFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);
3762: }
3763: } else {
3764: PetscCalloc2(Nf,&quads,Nf,&geoms);
3765: for (f = 0; f < Nf; ++f) {
3766: PetscObject obj;
3767: PetscClassId id;
3768: PetscBool fimp;
3770: PetscDSGetImplicit(prob, f, &fimp);
3771: if (isImplicit != fimp) continue;
3772: PetscDSGetDiscretization(prob, f, &obj);
3773: PetscObjectGetClassId(obj, &id);
3774: if (id == PETSCFE_CLASSID) {
3775: PetscFE fe = (PetscFE) obj;
3777: PetscFEGetQuadrature(fe, &quads[f]);
3778: PetscObjectReference((PetscObject)quads[f]);
3779: DMSNESGetFEGeom(coordField,cellIS,quads[f],PETSC_FALSE,&geoms[f]);
3780: }
3781: }
3782: }
3783: }
3784: /* Loop over chunks */
3785: ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
3786: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
3787: if (useFEM) {ISCreate(PETSC_COMM_SELF, &chunkIS);}
3788: numCells = cEnd - cStart;
3789: numChunks = 1;
3790: cellChunkSize = numCells/numChunks;
3791: numChunks = PetscMin(1,numCells);
3792: for (chunk = 0; chunk < numChunks; ++chunk) {
3793: PetscScalar *elemVec, *fluxL = NULL, *fluxR = NULL;
3794: PetscReal *vol = NULL;
3795: PetscFVFaceGeom *fgeom = NULL;
3796: PetscInt cS = cStart+chunk*cellChunkSize, cE = PetscMin(cS+cellChunkSize, cEnd), numCells = cE - cS, c;
3797: PetscInt numFaces = 0;
3799: /* Extract field coefficients */
3800: if (useFEM) {
3801: ISGetPointSubrange(chunkIS, cS, cE, cells);
3802: DMPlexGetCellFields(dm, chunkIS, locX, locX_t, locA, &u, &u_t, &a);
3803: DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);
3804: PetscArrayzero(elemVec, numCells*totDim);
3805: }
3806: /* TODO We will interlace both our field coefficients (u, u_t, uL, uR, etc.) and our output (elemVec, fL, fR). I think this works */
3807: /* Loop over fields */
3808: for (f = 0; f < Nf; ++f) {
3809: PetscObject obj;
3810: PetscClassId id;
3811: PetscBool fimp;
3812: PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;
3814: PetscDSGetImplicit(prob, f, &fimp);
3815: if (isImplicit != fimp) continue;
3816: PetscDSGetDiscretization(prob, f, &obj);
3817: PetscObjectGetClassId(obj, &id);
3818: if (id == PETSCFE_CLASSID) {
3819: PetscFE fe = (PetscFE) obj;
3820: PetscFEGeom *geom = affineGeom ? affineGeom : geoms[f];
3821: PetscFEGeom *chunkGeom = NULL;
3822: PetscQuadrature quad = affineQuad ? affineQuad : quads[f];
3823: PetscInt Nq, Nb;
3825: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
3826: PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
3827: PetscFEGetDimension(fe, &Nb);
3828: blockSize = Nb;
3829: batchSize = numBlocks * blockSize;
3830: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
3831: numChunks = numCells / (numBatches*batchSize);
3832: Ne = numChunks*numBatches*batchSize;
3833: Nr = numCells % (numBatches*batchSize);
3834: offset = numCells - Nr;
3835: /* Integrate FE residual to get elemVec (need fields at quadrature points) */
3836: /* 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) */
3837: PetscFEGeomGetChunk(geom,0,offset,&chunkGeom);
3838: PetscFEIntegrateResidual(prob, f, Ne, chunkGeom, u, u_t, probAux, a, t, elemVec);
3839: PetscFEGeomGetChunk(geom,offset,numCells,&chunkGeom);
3840: PetscFEIntegrateResidual(prob, f, Nr, chunkGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, &elemVec[offset*totDim]);
3841: PetscFEGeomRestoreChunk(geom,offset,numCells,&chunkGeom);
3842: } else if (id == PETSCFV_CLASSID) {
3843: PetscFV fv = (PetscFV) obj;
3845: Ne = numFaces;
3846: /* Riemann solve over faces (need fields at face centroids) */
3847: /* We need to evaluate FE fields at those coordinates */
3848: PetscFVIntegrateRHSFunction(fv, prob, f, Ne, fgeom, vol, uL, uR, fluxL, fluxR);
3849: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f);
3850: }
3851: /* Loop over domain */
3852: if (useFEM) {
3853: /* Add elemVec to locX */
3854: for (c = cS; c < cE; ++c) {
3855: const PetscInt cell = cells ? cells[c] : c;
3856: const PetscInt cind = c - cStart;
3858: if (mesh->printFEM > 1) {DMPrintCellVector(cell, name, totDim, &elemVec[cind*totDim]);}
3859: if (ghostLabel) {
3860: PetscInt ghostVal;
3862: DMLabelGetValue(ghostLabel,cell,&ghostVal);
3863: if (ghostVal > 0) continue;
3864: }
3865: DMPlexVecSetClosure(dm, section, locF, cell, &elemVec[cind*totDim], ADD_ALL_VALUES);
3866: }
3867: }
3868: /* Handle time derivative */
3869: if (locX_t) {
3870: PetscScalar *x_t, *fa;
3872: VecGetArray(locF, &fa);
3873: VecGetArray(locX_t, &x_t);
3874: for (f = 0; f < Nf; ++f) {
3875: PetscFV fv;
3876: PetscObject obj;
3877: PetscClassId id;
3878: PetscInt pdim, d;
3880: PetscDSGetDiscretization(prob, f, &obj);
3881: PetscObjectGetClassId(obj, &id);
3882: if (id != PETSCFV_CLASSID) continue;
3883: fv = (PetscFV) obj;
3884: PetscFVGetNumComponents(fv, &pdim);
3885: for (c = cS; c < cE; ++c) {
3886: const PetscInt cell = cells ? cells[c] : c;
3887: PetscScalar *u_t, *r;
3889: if (ghostLabel) {
3890: PetscInt ghostVal;
3892: DMLabelGetValue(ghostLabel, cell, &ghostVal);
3893: if (ghostVal > 0) continue;
3894: }
3895: DMPlexPointLocalFieldRead(dm, cell, f, x_t, &u_t);
3896: DMPlexPointLocalFieldRef(dm, cell, f, fa, &r);
3897: for (d = 0; d < pdim; ++d) r[d] += u_t[d];
3898: }
3899: }
3900: VecRestoreArray(locX_t, &x_t);
3901: VecRestoreArray(locF, &fa);
3902: }
3903: if (useFEM) {
3904: DMPlexRestoreCellFields(dm, chunkIS, locX, locX_t, locA, &u, &u_t, &a);
3905: DMRestoreWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);
3906: }
3907: }
3908: if (useFEM) {ISDestroy(&chunkIS);}
3909: ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);
3910: /* TODO Could include boundary residual here (see DMPlexComputeResidual_Internal) */
3911: if (useFEM) {
3912: if (maxDegree <= 1) {
3913: DMSNESRestoreFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);
3914: PetscQuadratureDestroy(&affineQuad);
3915: } else {
3916: for (f = 0; f < Nf; ++f) {
3917: DMSNESRestoreFEGeom(coordField,cellIS,quads[f],PETSC_FALSE,&geoms[f]);
3918: PetscQuadratureDestroy(&quads[f]);
3919: }
3920: PetscFree2(quads,geoms);
3921: }
3922: }
3923: PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);
3924: return(0);
3925: }
3927: /*
3928: 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
3930: X - The local solution vector
3931: X_t - The local solution time derviative vector, or NULL
3932: */
3933: PetscErrorCode DMPlexComputeJacobian_Patch_Internal(DM dm, PetscSection section, PetscSection globalSection, IS cellIS,
3934: PetscReal t, PetscReal X_tShift, Vec X, Vec X_t, Mat Jac, Mat JacP, void *ctx)
3935: {
3936: DM_Plex *mesh = (DM_Plex *) dm->data;
3937: const char *name = "Jacobian", *nameP = "JacobianPre";
3938: DM dmAux = NULL;
3939: PetscDS prob, probAux = NULL;
3940: PetscSection sectionAux = NULL;
3941: Vec A;
3942: DMField coordField;
3943: PetscFEGeom *cgeomFEM;
3944: PetscQuadrature qGeom = NULL;
3945: Mat J = Jac, JP = JacP;
3946: PetscScalar *work, *u = NULL, *u_t = NULL, *a = NULL, *elemMat = NULL, *elemMatP = NULL, *elemMatD = NULL;
3947: PetscBool hasJac, hasPrec, hasDyn, assembleJac, isMatIS, isMatISP, *isFE, hasFV = PETSC_FALSE;
3948: const PetscInt *cells;
3949: PetscInt Nf, fieldI, fieldJ, maxDegree, numCells, cStart, cEnd, numChunks, chunkSize, chunk, totDim, totDimAux = 0, sz, wsz, off = 0, offCell = 0;
3950: PetscErrorCode ierr;
3953: CHKMEMQ;
3954: ISGetLocalSize(cellIS, &numCells);
3955: ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
3956: PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);
3957: DMGetDS(dm, &prob);
3958: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
3959: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
3960: if (dmAux) {
3961: DMGetLocalSection(dmAux, §ionAux);
3962: DMGetDS(dmAux, &probAux);
3963: }
3964: /* Get flags */
3965: PetscDSGetNumFields(prob, &Nf);
3966: DMGetWorkArray(dm, Nf, MPIU_BOOL, &isFE);
3967: for (fieldI = 0; fieldI < Nf; ++fieldI) {
3968: PetscObject disc;
3969: PetscClassId id;
3970: PetscDSGetDiscretization(prob, fieldI, &disc);
3971: PetscObjectGetClassId(disc, &id);
3972: if (id == PETSCFE_CLASSID) {isFE[fieldI] = PETSC_TRUE;}
3973: else if (id == PETSCFV_CLASSID) {hasFV = PETSC_TRUE; isFE[fieldI] = PETSC_FALSE;}
3974: }
3975: PetscDSHasJacobian(prob, &hasJac);
3976: PetscDSHasJacobianPreconditioner(prob, &hasPrec);
3977: PetscDSHasDynamicJacobian(prob, &hasDyn);
3978: assembleJac = hasJac && hasPrec && (Jac != JacP) ? PETSC_TRUE : PETSC_FALSE;
3979: hasDyn = hasDyn && (X_tShift != 0.0) ? PETSC_TRUE : PETSC_FALSE;
3980: PetscObjectTypeCompare((PetscObject) Jac, MATIS, &isMatIS);
3981: PetscObjectTypeCompare((PetscObject) JacP, MATIS, &isMatISP);
3982: /* Setup input data and temp arrays (should be DMGetWorkArray) */
3983: if (isMatISP || isMatISP) {DMPlexGetSubdomainSection(dm, &globalSection);}
3984: if (isMatIS) {MatISGetLocalMat(Jac, &J);}
3985: if (isMatISP) {MatISGetLocalMat(JacP, &JP);}
3986: if (hasFV) {MatSetOption(JP, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);} /* No allocated space for FV stuff, so ignore the zero entries */
3987: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
3988: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
3989: PetscDSGetTotalDimension(prob, &totDim);
3990: if (probAux) {PetscDSGetTotalDimension(probAux, &totDimAux);}
3991: CHKMEMQ;
3992: /* Compute batch sizes */
3993: if (isFE[0]) {
3994: PetscFE fe;
3995: PetscQuadrature q;
3996: PetscInt numQuadPoints, numBatches, batchSize, numBlocks, blockSize, Nb;
3998: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
3999: PetscFEGetQuadrature(fe, &q);
4000: PetscQuadratureGetData(q, NULL, NULL, &numQuadPoints, NULL, NULL);
4001: PetscFEGetDimension(fe, &Nb);
4002: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
4003: blockSize = Nb*numQuadPoints;
4004: batchSize = numBlocks * blockSize;
4005: chunkSize = numBatches * batchSize;
4006: numChunks = numCells / chunkSize + numCells % chunkSize;
4007: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
4008: } else {
4009: chunkSize = numCells;
4010: numChunks = 1;
4011: }
4012: /* Get work space */
4013: wsz = (((X?1:0) + (X_t?1:0) + (dmAux?1:0))*totDim + ((hasJac?1:0) + (hasPrec?1:0) + (hasDyn?1:0))*totDim*totDim)*chunkSize;
4014: DMGetWorkArray(dm, wsz, MPIU_SCALAR, &work);
4015: PetscArrayzero(work, wsz);
4016: off = 0;
4017: u = X ? (sz = chunkSize*totDim, off += sz, work+off-sz) : NULL;
4018: u_t = X_t ? (sz = chunkSize*totDim, off += sz, work+off-sz) : NULL;
4019: a = dmAux ? (sz = chunkSize*totDimAux, off += sz, work+off-sz) : NULL;
4020: elemMat = hasJac ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL;
4021: elemMatP = hasPrec ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL;
4022: elemMatD = hasDyn ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL;
4023: if (off != wsz) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error is workspace size %D should be %D", off, wsz);
4024: /* Setup geometry */
4025: DMGetCoordinateField(dm, &coordField);
4026: DMFieldGetDegree(coordField, cellIS, NULL, &maxDegree);
4027: if (maxDegree <= 1) {DMFieldCreateDefaultQuadrature(coordField, cellIS, &qGeom);}
4028: if (!qGeom) {
4029: PetscFE fe;
4031: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
4032: PetscFEGetQuadrature(fe, &qGeom);
4033: PetscObjectReference((PetscObject) qGeom);
4034: }
4035: DMSNESGetFEGeom(coordField, cellIS, qGeom, PETSC_FALSE, &cgeomFEM);
4036: /* Compute volume integrals */
4037: if (assembleJac) {MatZeroEntries(J);}
4038: MatZeroEntries(JP);
4039: for (chunk = 0; chunk < numChunks; ++chunk, offCell += chunkSize) {
4040: const PetscInt Ncell = PetscMin(chunkSize, numCells - offCell);
4041: PetscInt c;
4043: /* Extract values */
4044: for (c = 0; c < Ncell; ++c) {
4045: const PetscInt cell = cells ? cells[c+offCell] : c+offCell;
4046: PetscScalar *x = NULL, *x_t = NULL;
4047: PetscInt i;
4049: if (X) {
4050: DMPlexVecGetClosure(dm, section, X, cell, NULL, &x);
4051: for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
4052: DMPlexVecRestoreClosure(dm, section, X, cell, NULL, &x);
4053: }
4054: if (X_t) {
4055: DMPlexVecGetClosure(dm, section, X_t, cell, NULL, &x_t);
4056: for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
4057: DMPlexVecRestoreClosure(dm, section, X_t, cell, NULL, &x_t);
4058: }
4059: if (dmAux) {
4060: DMPlexVecGetClosure(dmAux, sectionAux, A, cell, NULL, &x);
4061: for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
4062: DMPlexVecRestoreClosure(dmAux, sectionAux, A, cell, NULL, &x);
4063: }
4064: }
4065: CHKMEMQ;
4066: for (fieldI = 0; fieldI < Nf; ++fieldI) {
4067: PetscFE fe;
4068: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
4069: for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
4070: if (hasJac) {PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMat);}
4071: if (hasPrec) {PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_PRE, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMatP);}
4072: if (hasDyn) {PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMatD);}
4073: }
4074: /* For finite volume, add the identity */
4075: if (!isFE[fieldI]) {
4076: PetscFV fv;
4077: PetscInt eOffset = 0, Nc, fc, foff;
4079: PetscDSGetFieldOffset(prob, fieldI, &foff);
4080: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fv);
4081: PetscFVGetNumComponents(fv, &Nc);
4082: for (c = 0; c < chunkSize; ++c, eOffset += totDim*totDim) {
4083: for (fc = 0; fc < Nc; ++fc) {
4084: const PetscInt i = foff + fc;
4085: if (hasJac) {elemMat [eOffset+i*totDim+i] = 1.0;}
4086: if (hasPrec) {elemMatP[eOffset+i*totDim+i] = 1.0;}
4087: }
4088: }
4089: }
4090: }
4091: CHKMEMQ;
4092: /* Add contribution from X_t */
4093: if (hasDyn) {for (c = 0; c < chunkSize*totDim*totDim; ++c) elemMat[c] += X_tShift*elemMatD[c];}
4094: /* Insert values into matrix */
4095: for (c = 0; c < Ncell; ++c) {
4096: const PetscInt cell = cells ? cells[c+offCell] : c+offCell;
4097: if (mesh->printFEM > 1) {
4098: if (hasJac) {DMPrintCellMatrix(cell, name, totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);}
4099: if (hasPrec) {DMPrintCellMatrix(cell, nameP, totDim, totDim, &elemMatP[(c-cStart)*totDim*totDim]);}
4100: }
4101: if (assembleJac) {DMPlexMatSetClosure(dm, section, globalSection, Jac, cell, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);}
4102: DMPlexMatSetClosure(dm, section, globalSection, JP, cell, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);
4103: }
4104: CHKMEMQ;
4105: }
4106: /* Cleanup */
4107: DMSNESRestoreFEGeom(coordField, cellIS, qGeom, PETSC_FALSE, &cgeomFEM);
4108: PetscQuadratureDestroy(&qGeom);
4109: if (hasFV) {MatSetOption(JacP, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);}
4110: DMRestoreWorkArray(dm, Nf, MPIU_BOOL, &isFE);
4111: 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);
4112: /* Compute boundary integrals */
4113: /* DMPlexComputeBdJacobian_Internal(dm, X, X_t, t, X_tShift, Jac, JacP, ctx); */
4114: /* Assemble matrix */
4115: if (assembleJac) {MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);}
4116: MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);
4117: PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);
4118: CHKMEMQ;
4119: return(0);
4120: }
4122: /******** FEM Assembly Function ********/
4124: static PetscErrorCode DMConvertPlex_Internal(DM dm, DM *plex, PetscBool copy)
4125: {
4126: PetscBool isPlex;
4130: PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);
4131: if (isPlex) {
4132: *plex = dm;
4133: PetscObjectReference((PetscObject) dm);
4134: } else {
4135: PetscObjectQuery((PetscObject) dm, "dm_plex", (PetscObject *) plex);
4136: if (!*plex) {
4137: DMConvert(dm,DMPLEX,plex);
4138: PetscObjectCompose((PetscObject) dm, "dm_plex", (PetscObject) *plex);
4139: if (copy) {
4140: const char *comps[] = {"A", "dmAux"};
4141: PetscObject obj;
4142: PetscInt i;
4144: for (i = 0; i < 2; ++i) {
4145: PetscObjectQuery((PetscObject) dm, comps[i], &obj);
4146: PetscObjectCompose((PetscObject) *plex, comps[i], obj);
4147: }
4148: }
4149: } else {
4150: PetscObjectReference((PetscObject) *plex);
4151: }
4152: }
4153: return(0);
4154: }
4156: /*@
4157: DMPlexGetGeometryFVM - Return precomputed geometric data
4159: Collective on DM
4161: Input Parameter:
4162: . dm - The DM
4164: Output Parameters:
4165: + facegeom - The values precomputed from face geometry
4166: . cellgeom - The values precomputed from cell geometry
4167: - minRadius - The minimum radius over the mesh of an inscribed sphere in a cell
4169: Level: developer
4171: .seealso: DMPlexTSSetRHSFunctionLocal()
4172: @*/
4173: PetscErrorCode DMPlexGetGeometryFVM(DM dm, Vec *facegeom, Vec *cellgeom, PetscReal *minRadius)
4174: {
4175: DM plex;
4180: DMConvertPlex_Internal(dm,&plex,PETSC_TRUE);
4181: DMPlexGetDataFVM(plex, NULL, cellgeom, facegeom, NULL);
4182: if (minRadius) {DMPlexGetMinRadius(plex, minRadius);}
4183: DMDestroy(&plex);
4184: return(0);
4185: }
4187: /*@
4188: DMPlexGetGradientDM - Return gradient data layout
4190: Collective on DM
4192: Input Parameters:
4193: + dm - The DM
4194: - fv - The PetscFV
4196: Output Parameter:
4197: . dmGrad - The layout for gradient values
4199: Level: developer
4201: .seealso: DMPlexSNESGetGeometryFVM()
4202: @*/
4203: PetscErrorCode DMPlexGetGradientDM(DM dm, PetscFV fv, DM *dmGrad)
4204: {
4205: DM plex;
4206: PetscBool computeGradients;
4213: PetscFVGetComputeGradients(fv, &computeGradients);
4214: if (!computeGradients) {*dmGrad = NULL; return(0);}
4215: DMConvertPlex_Internal(dm,&plex,PETSC_TRUE);
4216: DMPlexGetDataFVM(plex, fv, NULL, NULL, dmGrad);
4217: DMDestroy(&plex);
4218: return(0);
4219: }
4221: static PetscErrorCode DMPlexComputeBdResidual_Single_Internal(DM dm, PetscReal t, DMLabel label, PetscInt numValues, const PetscInt values[], PetscInt field, Vec locX, Vec locX_t, Vec locF, DMField coordField, IS facetIS)
4222: {
4223: DM_Plex *mesh = (DM_Plex *) dm->data;
4224: DM plex = NULL, plexA = NULL;
4225: DMEnclosureType encAux;
4226: PetscDS prob, probAux = NULL;
4227: PetscSection section, sectionAux = NULL;
4228: Vec locA = NULL;
4229: PetscScalar *u = NULL, *u_t = NULL, *a = NULL, *elemVec = NULL;
4230: PetscInt v;
4231: PetscInt totDim, totDimAux = 0;
4232: PetscErrorCode ierr;
4235: DMConvert(dm, DMPLEX, &plex);
4236: DMGetLocalSection(dm, §ion);
4237: DMGetDS(dm, &prob);
4238: PetscDSGetTotalDimension(prob, &totDim);
4239: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
4240: if (locA) {
4241: DM dmAux;
4243: VecGetDM(locA, &dmAux);
4244: DMGetEnclosureRelation(dmAux, dm, &encAux);
4245: DMConvert(dmAux, DMPLEX, &plexA);
4246: DMGetDS(plexA, &probAux);
4247: PetscDSGetTotalDimension(probAux, &totDimAux);
4248: DMGetLocalSection(plexA, §ionAux);
4249: }
4250: for (v = 0; v < numValues; ++v) {
4251: PetscFEGeom *fgeom;
4252: PetscInt maxDegree;
4253: PetscQuadrature qGeom = NULL;
4254: IS pointIS;
4255: const PetscInt *points;
4256: PetscInt numFaces, face, Nq;
4258: DMLabelGetStratumIS(label, values[v], &pointIS);
4259: if (!pointIS) continue; /* No points with that id on this process */
4260: {
4261: IS isectIS;
4263: /* TODO: Special cases of ISIntersect where it is quick to check a priori if one is a superset of the other */
4264: ISIntersect_Caching_Internal(facetIS,pointIS,&isectIS);
4265: ISDestroy(&pointIS);
4266: pointIS = isectIS;
4267: }
4268: ISGetLocalSize(pointIS,&numFaces);
4269: ISGetIndices(pointIS,&points);
4270: PetscMalloc4(numFaces*totDim, &u, locX_t ? numFaces*totDim : 0, &u_t, numFaces*totDim, &elemVec, locA ? numFaces*totDimAux : 0, &a);
4271: DMFieldGetDegree(coordField,pointIS,NULL,&maxDegree);
4272: if (maxDegree <= 1) {
4273: DMFieldCreateDefaultQuadrature(coordField,pointIS,&qGeom);
4274: }
4275: if (!qGeom) {
4276: PetscFE fe;
4278: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
4279: PetscFEGetFaceQuadrature(fe, &qGeom);
4280: PetscObjectReference((PetscObject)qGeom);
4281: }
4282: PetscQuadratureGetData(qGeom, NULL, NULL, &Nq, NULL, NULL);
4283: DMSNESGetFEGeom(coordField,pointIS,qGeom,PETSC_TRUE,&fgeom);
4284: for (face = 0; face < numFaces; ++face) {
4285: const PetscInt point = points[face], *support;
4286: PetscScalar *x = NULL;
4287: PetscInt i;
4289: DMPlexGetSupport(dm, point, &support);
4290: DMPlexVecGetClosure(plex, section, locX, support[0], NULL, &x);
4291: for (i = 0; i < totDim; ++i) u[face*totDim+i] = x[i];
4292: DMPlexVecRestoreClosure(plex, section, locX, support[0], NULL, &x);
4293: if (locX_t) {
4294: DMPlexVecGetClosure(plex, section, locX_t, support[0], NULL, &x);
4295: for (i = 0; i < totDim; ++i) u_t[face*totDim+i] = x[i];
4296: DMPlexVecRestoreClosure(plex, section, locX_t, support[0], NULL, &x);
4297: }
4298: if (locA) {
4299: PetscInt subp;
4301: DMGetEnclosurePoint(plexA, dm, encAux, support[0], &subp);
4302: DMPlexVecGetClosure(plexA, sectionAux, locA, subp, NULL, &x);
4303: for (i = 0; i < totDimAux; ++i) a[face*totDimAux+i] = x[i];
4304: DMPlexVecRestoreClosure(plexA, sectionAux, locA, subp, NULL, &x);
4305: }
4306: }
4307: PetscArrayzero(elemVec, numFaces*totDim);
4308: {
4309: PetscFE fe;
4310: PetscInt Nb;
4311: PetscFEGeom *chunkGeom = NULL;
4312: /* Conforming batches */
4313: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
4314: /* Remainder */
4315: PetscInt Nr, offset;
4317: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
4318: PetscFEGetDimension(fe, &Nb);
4319: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
4320: /* TODO: documentation is unclear about what is going on with these numbers: how should Nb / Nq factor in ? */
4321: blockSize = Nb;
4322: batchSize = numBlocks * blockSize;
4323: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
4324: numChunks = numFaces / (numBatches*batchSize);
4325: Ne = numChunks*numBatches*batchSize;
4326: Nr = numFaces % (numBatches*batchSize);
4327: offset = numFaces - Nr;
4328: PetscFEGeomGetChunk(fgeom,0,offset,&chunkGeom);
4329: PetscFEIntegrateBdResidual(prob, field, Ne, chunkGeom, u, u_t, probAux, a, t, elemVec);
4330: PetscFEGeomRestoreChunk(fgeom, 0, offset, &chunkGeom);
4331: PetscFEGeomGetChunk(fgeom,offset,numFaces,&chunkGeom);
4332: PetscFEIntegrateBdResidual(prob, field, Nr, chunkGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, a ? &a[offset*totDimAux] : NULL, t, &elemVec[offset*totDim]);
4333: PetscFEGeomRestoreChunk(fgeom,offset,numFaces,&chunkGeom);
4334: }
4335: for (face = 0; face < numFaces; ++face) {
4336: const PetscInt point = points[face], *support;
4338: if (mesh->printFEM > 1) {DMPrintCellVector(point, "BdResidual", totDim, &elemVec[face*totDim]);}
4339: DMPlexGetSupport(plex, point, &support);
4340: DMPlexVecSetClosure(plex, NULL, locF, support[0], &elemVec[face*totDim], ADD_ALL_VALUES);
4341: }
4342: DMSNESRestoreFEGeom(coordField,pointIS,qGeom,PETSC_TRUE,&fgeom);
4343: PetscQuadratureDestroy(&qGeom);
4344: ISRestoreIndices(pointIS, &points);
4345: ISDestroy(&pointIS);
4346: PetscFree4(u, u_t, elemVec, a);
4347: }
4348: DMDestroy(&plex);
4349: DMDestroy(&plexA);
4350: return(0);
4351: }
4353: PetscErrorCode DMPlexComputeBdResidualSingle(DM dm, PetscReal t, DMLabel label, PetscInt numValues, const PetscInt values[], PetscInt field, Vec locX, Vec locX_t, Vec locF)
4354: {
4355: DMField coordField;
4356: DMLabel depthLabel;
4357: IS facetIS;
4358: PetscInt dim;
4362: DMGetDimension(dm, &dim);
4363: DMPlexGetDepthLabel(dm, &depthLabel);
4364: DMLabelGetStratumIS(depthLabel, dim-1, &facetIS);
4365: DMGetCoordinateField(dm, &coordField);
4366: DMPlexComputeBdResidual_Single_Internal(dm, t, label, numValues, values, field, locX, locX_t, locF, coordField, facetIS);
4367: ISDestroy(&facetIS);
4368: return(0);
4369: }
4371: PetscErrorCode DMPlexComputeBdResidual_Internal(DM dm, Vec locX, Vec locX_t, PetscReal t, Vec locF, void *user)
4372: {
4373: PetscDS prob;
4374: PetscInt numBd, bd;
4375: DMField coordField = NULL;
4376: IS facetIS = NULL;
4377: DMLabel depthLabel;
4378: PetscInt dim;
4382: DMGetDS(dm, &prob);
4383: DMPlexGetDepthLabel(dm, &depthLabel);
4384: DMGetDimension(dm, &dim);
4385: DMLabelGetStratumIS(depthLabel,dim - 1,&facetIS);
4386: PetscDSGetNumBoundary(prob, &numBd);
4387: for (bd = 0; bd < numBd; ++bd) {
4388: DMBoundaryConditionType type;
4389: const char *bdLabel;
4390: DMLabel label;
4391: const PetscInt *values;
4392: PetscInt field, numValues;
4393: PetscObject obj;
4394: PetscClassId id;
4396: PetscDSGetBoundary(prob, bd, &type, NULL, &bdLabel, &field, NULL, NULL, NULL, NULL, &numValues, &values, NULL);
4397: PetscDSGetDiscretization(prob, field, &obj);
4398: PetscObjectGetClassId(obj, &id);
4399: if ((id != PETSCFE_CLASSID) || (type & DM_BC_ESSENTIAL)) continue;
4400: if (!facetIS) {
4401: DMLabel depthLabel;
4402: PetscInt dim;
4404: DMPlexGetDepthLabel(dm, &depthLabel);
4405: DMGetDimension(dm, &dim);
4406: DMLabelGetStratumIS(depthLabel, dim - 1, &facetIS);
4407: }
4408: DMGetCoordinateField(dm, &coordField);
4409: DMGetLabel(dm, bdLabel, &label);
4410: DMPlexComputeBdResidual_Single_Internal(dm, t, label, numValues, values, field, locX, locX_t, locF, coordField, facetIS);
4411: }
4412: ISDestroy(&facetIS);
4413: return(0);
4414: }
4416: PetscErrorCode DMPlexComputeResidual_Internal(DM dm, IS cellIS, PetscReal time, Vec locX, Vec locX_t, PetscReal t, Vec locF, void *user)
4417: {
4418: DM_Plex *mesh = (DM_Plex *) dm->data;
4419: const char *name = "Residual";
4420: DM dmAux = NULL;
4421: DM dmGrad = NULL;
4422: DMLabel ghostLabel = NULL;
4423: PetscDS prob = NULL;
4424: PetscDS probAux = NULL;
4425: PetscSection section = NULL;
4426: PetscBool useFEM = PETSC_FALSE;
4427: PetscBool useFVM = PETSC_FALSE;
4428: PetscBool isImplicit = (locX_t || time == PETSC_MIN_REAL) ? PETSC_TRUE : PETSC_FALSE;
4429: PetscFV fvm = NULL;
4430: PetscFVCellGeom *cgeomFVM = NULL;
4431: PetscFVFaceGeom *fgeomFVM = NULL;
4432: DMField coordField = NULL;
4433: Vec locA, cellGeometryFVM = NULL, faceGeometryFVM = NULL, grad, locGrad = NULL;
4434: PetscScalar *u = NULL, *u_t, *a, *uL, *uR;
4435: IS chunkIS;
4436: const PetscInt *cells;
4437: PetscInt cStart, cEnd, numCells;
4438: PetscInt Nf, f, totDim, totDimAux, numChunks, cellChunkSize, faceChunkSize, chunk, fStart, fEnd;
4439: PetscInt maxDegree = PETSC_MAX_INT;
4440: PetscQuadrature affineQuad = NULL, *quads = NULL;
4441: PetscFEGeom *affineGeom = NULL, **geoms = NULL;
4442: PetscErrorCode ierr;
4445: PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);
4446: /* TODO The places where we have to use isFE are probably the member functions for the PetscDisc class */
4447: /* TODO The FVM geometry is over-manipulated. Make the precalc functions return exactly what we need */
4448: /* FEM+FVM */
4449: ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
4450: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
4451: /* 1: Get sizes from dm and dmAux */
4452: DMGetLocalSection(dm, §ion);
4453: DMGetLabel(dm, "ghost", &ghostLabel);
4454: DMGetCellDS(dm, cells ? cells[cStart] : cStart, &prob);
4455: PetscDSGetNumFields(prob, &Nf);
4456: PetscDSGetTotalDimension(prob, &totDim);
4457: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
4458: if (locA) {
4459: PetscInt subcell;
4460: VecGetDM(locA, &dmAux);
4461: DMGetEnclosurePoint(dmAux, dm, DM_ENC_UNKNOWN, cStart, &subcell);
4462: DMGetCellDS(dmAux, subcell, &probAux);
4463: PetscDSGetTotalDimension(probAux, &totDimAux);
4464: }
4465: /* 2: Get geometric data */
4466: for (f = 0; f < Nf; ++f) {
4467: PetscObject obj;
4468: PetscClassId id;
4469: PetscBool fimp;
4471: PetscDSGetImplicit(prob, f, &fimp);
4472: if (isImplicit != fimp) continue;
4473: PetscDSGetDiscretization(prob, f, &obj);
4474: PetscObjectGetClassId(obj, &id);
4475: if (id == PETSCFE_CLASSID) {useFEM = PETSC_TRUE;}
4476: if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;}
4477: }
4478: if (useFEM) {
4479: DMGetCoordinateField(dm, &coordField);
4480: DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);
4481: if (maxDegree <= 1) {
4482: DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);
4483: if (affineQuad) {
4484: DMSNESGetFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);
4485: }
4486: } else {
4487: PetscCalloc2(Nf,&quads,Nf,&geoms);
4488: for (f = 0; f < Nf; ++f) {
4489: PetscObject obj;
4490: PetscClassId id;
4491: PetscBool fimp;
4493: PetscDSGetImplicit(prob, f, &fimp);
4494: if (isImplicit != fimp) continue;
4495: PetscDSGetDiscretization(prob, f, &obj);
4496: PetscObjectGetClassId(obj, &id);
4497: if (id == PETSCFE_CLASSID) {
4498: PetscFE fe = (PetscFE) obj;
4500: PetscFEGetQuadrature(fe, &quads[f]);
4501: PetscObjectReference((PetscObject)quads[f]);
4502: DMSNESGetFEGeom(coordField,cellIS,quads[f],PETSC_FALSE,&geoms[f]);
4503: }
4504: }
4505: }
4506: }
4507: if (useFVM) {
4508: DMPlexGetGeometryFVM(dm, &faceGeometryFVM, &cellGeometryFVM, NULL);
4509: VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
4510: VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
4511: /* Reconstruct and limit cell gradients */
4512: DMPlexGetGradientDM(dm, fvm, &dmGrad);
4513: if (dmGrad) {
4514: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
4515: DMGetGlobalVector(dmGrad, &grad);
4516: DMPlexReconstructGradients_Internal(dm, fvm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
4517: /* Communicate gradient values */
4518: DMGetLocalVector(dmGrad, &locGrad);
4519: DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);
4520: DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);
4521: DMRestoreGlobalVector(dmGrad, &grad);
4522: }
4523: /* Handle non-essential (e.g. outflow) boundary values */
4524: DMPlexInsertBoundaryValues(dm, PETSC_FALSE, locX, time, faceGeometryFVM, cellGeometryFVM, locGrad);
4525: }
4526: /* Loop over chunks */
4527: if (useFEM) {ISCreate(PETSC_COMM_SELF, &chunkIS);}
4528: numCells = cEnd - cStart;
4529: numChunks = 1;
4530: cellChunkSize = numCells/numChunks;
4531: faceChunkSize = (fEnd - fStart)/numChunks;
4532: numChunks = PetscMin(1,numCells);
4533: for (chunk = 0; chunk < numChunks; ++chunk) {
4534: PetscScalar *elemVec, *fluxL, *fluxR;
4535: PetscReal *vol;
4536: PetscFVFaceGeom *fgeom;
4537: PetscInt cS = cStart+chunk*cellChunkSize, cE = PetscMin(cS+cellChunkSize, cEnd), numCells = cE - cS, c;
4538: PetscInt fS = fStart+chunk*faceChunkSize, fE = PetscMin(fS+faceChunkSize, fEnd), numFaces = 0, face;
4540: /* Extract field coefficients */
4541: if (useFEM) {
4542: ISGetPointSubrange(chunkIS, cS, cE, cells);
4543: DMPlexGetCellFields(dm, chunkIS, locX, locX_t, locA, &u, &u_t, &a);
4544: DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);
4545: PetscArrayzero(elemVec, numCells*totDim);
4546: }
4547: if (useFVM) {
4548: DMPlexGetFaceFields(dm, fS, fE, locX, locX_t, faceGeometryFVM, cellGeometryFVM, locGrad, &numFaces, &uL, &uR);
4549: DMPlexGetFaceGeometry(dm, fS, fE, faceGeometryFVM, cellGeometryFVM, &numFaces, &fgeom, &vol);
4550: DMGetWorkArray(dm, numFaces*totDim, MPIU_SCALAR, &fluxL);
4551: DMGetWorkArray(dm, numFaces*totDim, MPIU_SCALAR, &fluxR);
4552: PetscArrayzero(fluxL, numFaces*totDim);
4553: PetscArrayzero(fluxR, numFaces*totDim);
4554: }
4555: /* TODO We will interlace both our field coefficients (u, u_t, uL, uR, etc.) and our output (elemVec, fL, fR). I think this works */
4556: /* Loop over fields */
4557: for (f = 0; f < Nf; ++f) {
4558: PetscObject obj;
4559: PetscClassId id;
4560: PetscBool fimp;
4561: PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;
4563: PetscDSGetImplicit(prob, f, &fimp);
4564: if (isImplicit != fimp) continue;
4565: PetscDSGetDiscretization(prob, f, &obj);
4566: PetscObjectGetClassId(obj, &id);
4567: if (id == PETSCFE_CLASSID) {
4568: PetscFE fe = (PetscFE) obj;
4569: PetscFEGeom *geom = affineGeom ? affineGeom : geoms[f];
4570: PetscFEGeom *chunkGeom = NULL;
4571: PetscQuadrature quad = affineQuad ? affineQuad : quads[f];
4572: PetscInt Nq, Nb;
4574: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
4575: PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
4576: PetscFEGetDimension(fe, &Nb);
4577: blockSize = Nb;
4578: batchSize = numBlocks * blockSize;
4579: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
4580: numChunks = numCells / (numBatches*batchSize);
4581: Ne = numChunks*numBatches*batchSize;
4582: Nr = numCells % (numBatches*batchSize);
4583: offset = numCells - Nr;
4584: /* Integrate FE residual to get elemVec (need fields at quadrature points) */
4585: /* 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) */
4586: PetscFEGeomGetChunk(geom,0,offset,&chunkGeom);
4587: PetscFEIntegrateResidual(prob, f, Ne, chunkGeom, u, u_t, probAux, a, t, elemVec);
4588: PetscFEGeomGetChunk(geom,offset,numCells,&chunkGeom);
4589: PetscFEIntegrateResidual(prob, f, Nr, chunkGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, &elemVec[offset*totDim]);
4590: PetscFEGeomRestoreChunk(geom,offset,numCells,&chunkGeom);
4591: } else if (id == PETSCFV_CLASSID) {
4592: PetscFV fv = (PetscFV) obj;
4594: Ne = numFaces;
4595: /* Riemann solve over faces (need fields at face centroids) */
4596: /* We need to evaluate FE fields at those coordinates */
4597: PetscFVIntegrateRHSFunction(fv, prob, f, Ne, fgeom, vol, uL, uR, fluxL, fluxR);
4598: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f);
4599: }
4600: /* Loop over domain */
4601: if (useFEM) {
4602: /* Add elemVec to locX */
4603: for (c = cS; c < cE; ++c) {
4604: const PetscInt cell = cells ? cells[c] : c;
4605: const PetscInt cind = c - cStart;
4607: if (mesh->printFEM > 1) {DMPrintCellVector(cell, name, totDim, &elemVec[cind*totDim]);}
4608: if (ghostLabel) {
4609: PetscInt ghostVal;
4611: DMLabelGetValue(ghostLabel,cell,&ghostVal);
4612: if (ghostVal > 0) continue;
4613: }
4614: DMPlexVecSetClosure(dm, section, locF, cell, &elemVec[cind*totDim], ADD_ALL_VALUES);
4615: }
4616: }
4617: if (useFVM) {
4618: PetscScalar *fa;
4619: PetscInt iface;
4621: VecGetArray(locF, &fa);
4622: for (f = 0; f < Nf; ++f) {
4623: PetscFV fv;
4624: PetscObject obj;
4625: PetscClassId id;
4626: PetscInt foff, pdim;
4628: PetscDSGetDiscretization(prob, f, &obj);
4629: PetscDSGetFieldOffset(prob, f, &foff);
4630: PetscObjectGetClassId(obj, &id);
4631: if (id != PETSCFV_CLASSID) continue;
4632: fv = (PetscFV) obj;
4633: PetscFVGetNumComponents(fv, &pdim);
4634: /* Accumulate fluxes to cells */
4635: for (face = fS, iface = 0; face < fE; ++face) {
4636: const PetscInt *scells;
4637: PetscScalar *fL = NULL, *fR = NULL;
4638: PetscInt ghost, d, nsupp, nchild;
4640: DMLabelGetValue(ghostLabel, face, &ghost);
4641: DMPlexGetSupportSize(dm, face, &nsupp);
4642: DMPlexGetTreeChildren(dm, face, &nchild, NULL);
4643: if (ghost >= 0 || nsupp > 2 || nchild > 0) continue;
4644: DMPlexGetSupport(dm, face, &scells);
4645: DMLabelGetValue(ghostLabel,scells[0],&ghost);
4646: if (ghost <= 0) {DMPlexPointLocalFieldRef(dm, scells[0], f, fa, &fL);}
4647: DMLabelGetValue(ghostLabel,scells[1],&ghost);
4648: if (ghost <= 0) {DMPlexPointLocalFieldRef(dm, scells[1], f, fa, &fR);}
4649: for (d = 0; d < pdim; ++d) {
4650: if (fL) fL[d] -= fluxL[iface*totDim+foff+d];
4651: if (fR) fR[d] += fluxR[iface*totDim+foff+d];
4652: }
4653: ++iface;
4654: }
4655: }
4656: VecRestoreArray(locF, &fa);
4657: }
4658: /* Handle time derivative */
4659: if (locX_t) {
4660: PetscScalar *x_t, *fa;
4662: VecGetArray(locF, &fa);
4663: VecGetArray(locX_t, &x_t);
4664: for (f = 0; f < Nf; ++f) {
4665: PetscFV fv;
4666: PetscObject obj;
4667: PetscClassId id;
4668: PetscInt pdim, d;
4670: PetscDSGetDiscretization(prob, f, &obj);
4671: PetscObjectGetClassId(obj, &id);
4672: if (id != PETSCFV_CLASSID) continue;
4673: fv = (PetscFV) obj;
4674: PetscFVGetNumComponents(fv, &pdim);
4675: for (c = cS; c < cE; ++c) {
4676: const PetscInt cell = cells ? cells[c] : c;
4677: PetscScalar *u_t, *r;
4679: if (ghostLabel) {
4680: PetscInt ghostVal;
4682: DMLabelGetValue(ghostLabel, cell, &ghostVal);
4683: if (ghostVal > 0) continue;
4684: }
4685: DMPlexPointLocalFieldRead(dm, cell, f, x_t, &u_t);
4686: DMPlexPointLocalFieldRef(dm, cell, f, fa, &r);
4687: for (d = 0; d < pdim; ++d) r[d] += u_t[d];
4688: }
4689: }
4690: VecRestoreArray(locX_t, &x_t);
4691: VecRestoreArray(locF, &fa);
4692: }
4693: if (useFEM) {
4694: DMPlexRestoreCellFields(dm, chunkIS, locX, locX_t, locA, &u, &u_t, &a);
4695: DMRestoreWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);
4696: }
4697: if (useFVM) {
4698: DMPlexRestoreFaceFields(dm, fS, fE, locX, locX_t, faceGeometryFVM, cellGeometryFVM, locGrad, &numFaces, &uL, &uR);
4699: DMPlexRestoreFaceGeometry(dm, fS, fE, faceGeometryFVM, cellGeometryFVM, &numFaces, &fgeom, &vol);
4700: DMRestoreWorkArray(dm, numFaces*totDim, MPIU_SCALAR, &fluxL);
4701: DMRestoreWorkArray(dm, numFaces*totDim, MPIU_SCALAR, &fluxR);
4702: if (dmGrad) {DMRestoreLocalVector(dmGrad, &locGrad);}
4703: }
4704: }
4705: if (useFEM) {ISDestroy(&chunkIS);}
4706: ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);
4708: if (useFEM) {
4709: DMPlexComputeBdResidual_Internal(dm, locX, locX_t, t, locF, user);
4711: if (maxDegree <= 1) {
4712: DMSNESRestoreFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);
4713: PetscQuadratureDestroy(&affineQuad);
4714: } else {
4715: for (f = 0; f < Nf; ++f) {
4716: DMSNESRestoreFEGeom(coordField,cellIS,quads[f],PETSC_FALSE,&geoms[f]);
4717: PetscQuadratureDestroy(&quads[f]);
4718: }
4719: PetscFree2(quads,geoms);
4720: }
4721: }
4723: /* FEM */
4724: /* 1: Get sizes from dm and dmAux */
4725: /* 2: Get geometric data */
4726: /* 3: Handle boundary values */
4727: /* 4: Loop over domain */
4728: /* Extract coefficients */
4729: /* Loop over fields */
4730: /* Set tiling for FE*/
4731: /* Integrate FE residual to get elemVec */
4732: /* Loop over subdomain */
4733: /* Loop over quad points */
4734: /* Transform coords to real space */
4735: /* Evaluate field and aux fields at point */
4736: /* Evaluate residual at point */
4737: /* Transform residual to real space */
4738: /* Add residual to elemVec */
4739: /* Loop over domain */
4740: /* Add elemVec to locX */
4742: /* FVM */
4743: /* Get geometric data */
4744: /* If using gradients */
4745: /* Compute gradient data */
4746: /* Loop over domain faces */
4747: /* Count computational faces */
4748: /* Reconstruct cell gradient */
4749: /* Loop over domain cells */
4750: /* Limit cell gradients */
4751: /* Handle boundary values */
4752: /* Loop over domain faces */
4753: /* Read out field, centroid, normal, volume for each side of face */
4754: /* Riemann solve over faces */
4755: /* Loop over domain faces */
4756: /* Accumulate fluxes to cells */
4757: /* TODO Change printFEM to printDisc here */
4758: if (mesh->printFEM) {
4759: Vec locFbc;
4760: PetscInt pStart, pEnd, p, maxDof;
4761: PetscScalar *zeroes;
4763: VecDuplicate(locF,&locFbc);
4764: VecCopy(locF,locFbc);
4765: PetscSectionGetChart(section,&pStart,&pEnd);
4766: PetscSectionGetMaxDof(section,&maxDof);
4767: PetscCalloc1(maxDof,&zeroes);
4768: for (p = pStart; p < pEnd; p++) {
4769: VecSetValuesSection(locFbc,section,p,zeroes,INSERT_BC_VALUES);
4770: }
4771: PetscFree(zeroes);
4772: DMPrintLocalVec(dm, name, mesh->printTol, locFbc);
4773: VecDestroy(&locFbc);
4774: }
4775: PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);
4776: return(0);
4777: }
4779: PetscErrorCode DMPlexComputeResidual_Hybrid_Internal(DM dm, IS cellIS, PetscReal time, Vec locX, Vec locX_t, PetscReal t, Vec locF, void *user)
4780: {
4781: DM_Plex *mesh = (DM_Plex *) dm->data;
4782: const char *name = "Hybrid Residual";
4783: DM dmAux = NULL;
4784: DMLabel ghostLabel = NULL;
4785: PetscDS prob = NULL;
4786: PetscDS probAux = NULL;
4787: PetscSection section = NULL;
4788: DMField coordField = NULL;
4789: Vec locA;
4790: PetscScalar *u = NULL, *u_t, *a;
4791: PetscScalar *elemVec;
4792: IS chunkIS;
4793: const PetscInt *cells;
4794: PetscInt *faces;
4795: PetscInt cStart, cEnd, numCells;
4796: PetscInt Nf, f, totDim, totDimAux, numChunks, cellChunkSize, chunk;
4797: PetscInt maxDegree = PETSC_MAX_INT;
4798: PetscQuadrature affineQuad = NULL, *quads = NULL;
4799: PetscFEGeom *affineGeom = NULL, **geoms = NULL;
4800: PetscErrorCode ierr;
4803: PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);
4804: /* TODO The places where we have to use isFE are probably the member functions for the PetscDisc class */
4805: /* FEM */
4806: ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
4807: /* 1: Get sizes from dm and dmAux */
4808: DMGetSection(dm, §ion);
4809: DMGetLabel(dm, "ghost", &ghostLabel);
4810: DMGetCellDS(dm, cStart, &prob);
4811: PetscDSGetNumFields(prob, &Nf);
4812: PetscDSGetTotalDimension(prob, &totDim);
4813: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
4814: if (locA) {
4815: VecGetDM(locA, &dmAux);
4816: DMGetCellDS(dmAux, cStart, &probAux);
4817: PetscDSGetTotalDimension(probAux, &totDimAux);
4818: }
4819: /* 2: Setup geometric data */
4820: DMGetCoordinateField(dm, &coordField);
4821: DMFieldGetDegree(coordField, cellIS, NULL, &maxDegree);
4822: if (maxDegree > 1) {
4823: PetscCalloc2(Nf,&quads,Nf,&geoms);
4824: for (f = 0; f < Nf; ++f) {
4825: PetscFE fe;
4827: PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
4828: if (fe) {
4829: PetscFEGetQuadrature(fe, &quads[f]);
4830: PetscObjectReference((PetscObject) quads[f]);
4831: }
4832: }
4833: }
4834: /* Loop over chunks */
4835: numCells = cEnd - cStart;
4836: cellChunkSize = numCells;
4837: numChunks = !numCells ? 0 : PetscCeilReal(((PetscReal) numCells)/cellChunkSize);
4838: PetscCalloc1(2*cellChunkSize, &faces);
4839: ISCreateGeneral(PETSC_COMM_SELF, cellChunkSize, faces, PETSC_USE_POINTER, &chunkIS);
4840: /* Extract field coefficients */
4841: /* NOTE This needs the end cap faces to have identical orientations */
4842: DMPlexGetCellFields(dm, cellIS, locX, locX_t, locA, &u, &u_t, &a);
4843: DMGetWorkArray(dm, cellChunkSize*totDim, MPIU_SCALAR, &elemVec);
4844: for (chunk = 0; chunk < numChunks; ++chunk) {
4845: PetscInt cS = cStart+chunk*cellChunkSize, cE = PetscMin(cS+cellChunkSize, cEnd), numCells = cE - cS, c;
4847: PetscMemzero(elemVec, cellChunkSize*totDim * sizeof(PetscScalar));
4848: /* Get faces */
4849: for (c = cS; c < cE; ++c) {
4850: const PetscInt cell = cells ? cells[c] : c;
4851: const PetscInt *cone;
4852: DMPlexGetCone(dm, cell, &cone);
4853: faces[(c-cS)*2+0] = cone[0];
4854: faces[(c-cS)*2+1] = cone[1];
4855: }
4856: ISGeneralSetIndices(chunkIS, cellChunkSize, faces, PETSC_USE_POINTER);
4857: /* Get geometric data */
4858: if (maxDegree <= 1) {
4859: if (!affineQuad) {DMFieldCreateDefaultQuadrature(coordField, chunkIS, &affineQuad);}
4860: if (affineQuad) {DMSNESGetFEGeom(coordField, chunkIS, affineQuad, PETSC_TRUE, &affineGeom);}
4861: } else {
4862: for (f = 0; f < Nf; ++f) {
4863: if (quads[f]) {DMSNESGetFEGeom(coordField, chunkIS, quads[f], PETSC_TRUE, &geoms[f]);}
4864: }
4865: }
4866: /* Loop over fields */
4867: for (f = 0; f < Nf; ++f) {
4868: PetscFE fe;
4869: PetscFEGeom *geom = affineGeom ? affineGeom : geoms[f];
4870: PetscFEGeom *chunkGeom = NULL;
4871: PetscQuadrature quad = affineQuad ? affineQuad : quads[f];
4872: PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset, Nq, Nb;
4874: PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
4875: if (!fe) continue;
4876: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
4877: PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
4878: PetscFEGetDimension(fe, &Nb);
4879: blockSize = Nb;
4880: batchSize = numBlocks * blockSize;
4881: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
4882: numChunks = numCells / (numBatches*batchSize);
4883: Ne = numChunks*numBatches*batchSize;
4884: Nr = numCells % (numBatches*batchSize);
4885: offset = numCells - Nr;
4886: PetscFEGeomGetChunk(geom,0,offset,&chunkGeom);
4887: PetscFEIntegrateHybridResidual(prob, f, Ne, chunkGeom, u, u_t, probAux, a, t, elemVec);
4888: PetscFEGeomGetChunk(geom,offset,numCells,&chunkGeom);
4889: PetscFEIntegrateHybridResidual(prob, f, Nr, chunkGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, &elemVec[offset*totDim]);
4890: PetscFEGeomRestoreChunk(geom,offset,numCells,&chunkGeom);
4891: }
4892: /* Add elemVec to locX */
4893: for (c = cS; c < cE; ++c) {
4894: const PetscInt cell = cells ? cells[c] : c;
4895: const PetscInt cind = c - cStart;
4897: if (mesh->printFEM > 1) {DMPrintCellVector(cell, name, totDim, &elemVec[cind*totDim]);}
4898: if (ghostLabel) {
4899: PetscInt ghostVal;
4901: DMLabelGetValue(ghostLabel,cell,&ghostVal);
4902: if (ghostVal > 0) continue;
4903: }
4904: DMPlexVecSetClosure(dm, section, locF, cell, &elemVec[cind*totDim], ADD_ALL_VALUES);
4905: }
4906: }
4907: DMPlexRestoreCellFields(dm, cellIS, locX, locX_t, locA, &u, &u_t, &a);
4908: DMRestoreWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);
4909: PetscFree(faces);
4910: ISDestroy(&chunkIS);
4911: ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);
4912: if (maxDegree <= 1) {
4913: DMSNESRestoreFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);
4914: PetscQuadratureDestroy(&affineQuad);
4915: } else {
4916: for (f = 0; f < Nf; ++f) {
4917: if (geoms) {DMSNESRestoreFEGeom(coordField,cellIS,quads[f],PETSC_FALSE,&geoms[f]);}
4918: if (quads) {PetscQuadratureDestroy(&quads[f]);}
4919: }
4920: PetscFree2(quads,geoms);
4921: }
4922: PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);
4923: return(0);
4924: }
4926: PetscErrorCode DMPlexComputeBdJacobian_Single_Internal(DM dm, PetscReal t, DMLabel label, PetscInt numValues, const PetscInt values[], PetscInt fieldI, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, DMField coordField, IS facetIS)
4927: {
4928: DM_Plex *mesh = (DM_Plex *) dm->data;
4929: DM plex = NULL, plexA = NULL, tdm;
4930: DMEnclosureType encAux;
4931: PetscDS prob, probAux = NULL;
4932: PetscSection section, sectionAux = NULL;
4933: PetscSection globalSection, subSection = NULL;
4934: Vec locA = NULL, tv;
4935: PetscScalar *u = NULL, *u_t = NULL, *a = NULL, *elemMat = NULL;
4936: PetscInt v;
4937: PetscInt Nf, totDim, totDimAux = 0;
4938: PetscBool isMatISP, transform;
4939: PetscErrorCode ierr;
4942: DMConvert(dm, DMPLEX, &plex);
4943: DMHasBasisTransform(dm, &transform);
4944: DMGetBasisTransformDM_Internal(dm, &tdm);
4945: DMGetBasisTransformVec_Internal(dm, &tv);
4946: DMGetLocalSection(dm, §ion);
4947: DMGetDS(dm, &prob);
4948: PetscDSGetNumFields(prob, &Nf);
4949: PetscDSGetTotalDimension(prob, &totDim);
4950: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
4951: if (locA) {
4952: DM dmAux;
4954: VecGetDM(locA, &dmAux);
4955: DMGetEnclosureRelation(dmAux, dm, &encAux);
4956: DMConvert(dmAux, DMPLEX, &plexA);
4957: DMGetDS(plexA, &probAux);
4958: PetscDSGetTotalDimension(probAux, &totDimAux);
4959: DMGetLocalSection(plexA, §ionAux);
4960: }
4962: PetscObjectTypeCompare((PetscObject) JacP, MATIS, &isMatISP);
4963: DMGetGlobalSection(dm, &globalSection);
4964: if (isMatISP) {DMPlexGetSubdomainSection(dm, &subSection);}
4965: for (v = 0; v < numValues; ++v) {
4966: PetscFEGeom *fgeom;
4967: PetscInt maxDegree;
4968: PetscQuadrature qGeom = NULL;
4969: IS pointIS;
4970: const PetscInt *points;
4971: PetscInt numFaces, face, Nq;
4973: DMLabelGetStratumIS(label, values[v], &pointIS);
4974: if (!pointIS) continue; /* No points with that id on this process */
4975: {
4976: IS isectIS;
4978: /* TODO: Special cases of ISIntersect where it is quick to check a prior if one is a superset of the other */
4979: ISIntersect_Caching_Internal(facetIS,pointIS,&isectIS);
4980: ISDestroy(&pointIS);
4981: pointIS = isectIS;
4982: }
4983: ISGetLocalSize(pointIS, &numFaces);
4984: ISGetIndices(pointIS, &points);
4985: PetscMalloc4(numFaces*totDim, &u, locX_t ? numFaces*totDim : 0, &u_t, numFaces*totDim*totDim, &elemMat, locA ? numFaces*totDimAux : 0, &a);
4986: DMFieldGetDegree(coordField,pointIS,NULL,&maxDegree);
4987: if (maxDegree <= 1) {
4988: DMFieldCreateDefaultQuadrature(coordField,pointIS,&qGeom);
4989: }
4990: if (!qGeom) {
4991: PetscFE fe;
4993: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
4994: PetscFEGetFaceQuadrature(fe, &qGeom);
4995: PetscObjectReference((PetscObject)qGeom);
4996: }
4997: PetscQuadratureGetData(qGeom, NULL, NULL, &Nq, NULL, NULL);
4998: DMSNESGetFEGeom(coordField,pointIS,qGeom,PETSC_TRUE,&fgeom);
4999: for (face = 0; face < numFaces; ++face) {
5000: const PetscInt point = points[face], *support;
5001: PetscScalar *x = NULL;
5002: PetscInt i;
5004: DMPlexGetSupport(dm, point, &support);
5005: DMPlexVecGetClosure(plex, section, locX, support[0], NULL, &x);
5006: for (i = 0; i < totDim; ++i) u[face*totDim+i] = x[i];
5007: DMPlexVecRestoreClosure(plex, section, locX, support[0], NULL, &x);
5008: if (locX_t) {
5009: DMPlexVecGetClosure(plex, section, locX_t, support[0], NULL, &x);
5010: for (i = 0; i < totDim; ++i) u_t[face*totDim+i] = x[i];
5011: DMPlexVecRestoreClosure(plex, section, locX_t, support[0], NULL, &x);
5012: }
5013: if (locA) {
5014: PetscInt subp;
5015: DMGetEnclosurePoint(plexA, dm, encAux, support[0], &subp);
5016: DMPlexVecGetClosure(plexA, sectionAux, locA, subp, NULL, &x);
5017: for (i = 0; i < totDimAux; ++i) a[face*totDimAux+i] = x[i];
5018: DMPlexVecRestoreClosure(plexA, sectionAux, locA, subp, NULL, &x);
5019: }
5020: }
5021: PetscArrayzero(elemMat, numFaces*totDim*totDim);
5022: {
5023: PetscFE fe;
5024: PetscInt Nb;
5025: /* Conforming batches */
5026: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
5027: /* Remainder */
5028: PetscFEGeom *chunkGeom = NULL;
5029: PetscInt fieldJ, Nr, offset;
5031: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
5032: PetscFEGetDimension(fe, &Nb);
5033: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
5034: blockSize = Nb;
5035: batchSize = numBlocks * blockSize;
5036: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
5037: numChunks = numFaces / (numBatches*batchSize);
5038: Ne = numChunks*numBatches*batchSize;
5039: Nr = numFaces % (numBatches*batchSize);
5040: offset = numFaces - Nr;
5041: PetscFEGeomGetChunk(fgeom,0,offset,&chunkGeom);
5042: for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
5043: PetscFEIntegrateBdJacobian(prob, fieldI, fieldJ, Ne, chunkGeom, u, u_t, probAux, a, t, X_tShift, elemMat);
5044: }
5045: PetscFEGeomGetChunk(fgeom,offset,numFaces,&chunkGeom);
5046: for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
5047: PetscFEIntegrateBdJacobian(prob, fieldI, fieldJ, Nr, chunkGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, a ? &a[offset*totDimAux] : NULL, t, X_tShift, &elemMat[offset*totDim*totDim]);
5048: }
5049: PetscFEGeomRestoreChunk(fgeom,offset,numFaces,&chunkGeom);
5050: }
5051: for (face = 0; face < numFaces; ++face) {
5052: const PetscInt point = points[face], *support;
5054: /* Transform to global basis before insertion in Jacobian */
5055: DMPlexGetSupport(plex, point, &support);
5056: if (transform) {DMPlexBasisTransformPointTensor_Internal(dm, tdm, tv, support[0], PETSC_TRUE, totDim, &elemMat[face*totDim*totDim]);}
5057: if (mesh->printFEM > 1) {DMPrintCellMatrix(point, "BdJacobian", totDim, totDim, &elemMat[face*totDim*totDim]);}
5058: if (!isMatISP) {
5059: DMPlexMatSetClosure(plex, section, globalSection, JacP, support[0], &elemMat[face*totDim*totDim], ADD_VALUES);
5060: } else {
5061: Mat lJ;
5063: MatISGetLocalMat(JacP, &lJ);
5064: DMPlexMatSetClosure(plex, section, subSection, lJ, support[0], &elemMat[face*totDim*totDim], ADD_VALUES);
5065: }
5066: }
5067: DMSNESRestoreFEGeom(coordField,pointIS,qGeom,PETSC_TRUE,&fgeom);
5068: PetscQuadratureDestroy(&qGeom);
5069: ISRestoreIndices(pointIS, &points);
5070: ISDestroy(&pointIS);
5071: PetscFree4(u, u_t, elemMat, a);
5072: }
5073: if (plex) {DMDestroy(&plex);}
5074: if (plexA) {DMDestroy(&plexA);}
5075: return(0);
5076: }
5078: PetscErrorCode DMPlexComputeBdJacobianSingle(DM dm, PetscReal t, DMLabel label, PetscInt numValues, const PetscInt values[], PetscInt field, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP)
5079: {
5080: DMField coordField;
5081: DMLabel depthLabel;
5082: IS facetIS;
5083: PetscInt dim;
5087: DMGetDimension(dm, &dim);
5088: DMPlexGetDepthLabel(dm, &depthLabel);
5089: DMLabelGetStratumIS(depthLabel, dim-1, &facetIS);
5090: DMGetCoordinateField(dm, &coordField);
5091: DMPlexComputeBdJacobian_Single_Internal(dm, t, label, numValues, values, field, locX, locX_t, X_tShift, Jac, JacP, coordField, facetIS);
5092: ISDestroy(&facetIS);
5093: return(0);
5094: }
5096: PetscErrorCode DMPlexComputeBdJacobian_Internal(DM dm, Vec locX, Vec locX_t, PetscReal t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user)
5097: {
5098: PetscDS prob;
5099: PetscInt dim, numBd, bd;
5100: DMLabel depthLabel;
5101: DMField coordField = NULL;
5102: IS facetIS;
5103: PetscErrorCode ierr;
5106: DMGetDS(dm, &prob);
5107: DMPlexGetDepthLabel(dm, &depthLabel);
5108: DMGetDimension(dm, &dim);
5109: DMLabelGetStratumIS(depthLabel, dim-1, &facetIS);
5110: PetscDSGetNumBoundary(prob, &numBd);
5111: DMGetCoordinateField(dm, &coordField);
5112: for (bd = 0; bd < numBd; ++bd) {
5113: DMBoundaryConditionType type;
5114: const char *bdLabel;
5115: DMLabel label;
5116: const PetscInt *values;
5117: PetscInt fieldI, numValues;
5118: PetscObject obj;
5119: PetscClassId id;
5121: PetscDSGetBoundary(prob, bd, &type, NULL, &bdLabel, &fieldI, NULL, NULL, NULL, NULL, &numValues, &values, NULL);
5122: PetscDSGetDiscretization(prob, fieldI, &obj);
5123: PetscObjectGetClassId(obj, &id);
5124: if ((id != PETSCFE_CLASSID) || (type & DM_BC_ESSENTIAL)) continue;
5125: DMGetLabel(dm, bdLabel, &label);
5126: DMPlexComputeBdJacobian_Single_Internal(dm, t, label, numValues, values, fieldI, locX, locX_t, X_tShift, Jac, JacP, coordField, facetIS);
5127: }
5128: ISDestroy(&facetIS);
5129: return(0);
5130: }
5132: PetscErrorCode DMPlexComputeJacobian_Internal(DM dm, IS cellIS, PetscReal t, PetscReal X_tShift, Vec X, Vec X_t, Mat Jac, Mat JacP,void *user)
5133: {
5134: DM_Plex *mesh = (DM_Plex *) dm->data;
5135: const char *name = "Jacobian";
5136: DM dmAux, plex, tdm;
5137: DMEnclosureType encAux;
5138: Vec A, tv;
5139: DMField coordField;
5140: PetscDS prob, probAux = NULL;
5141: PetscSection section, globalSection, subSection, sectionAux;
5142: PetscScalar *elemMat, *elemMatP, *elemMatD, *u, *u_t, *a = NULL;
5143: const PetscInt *cells;
5144: PetscInt Nf, fieldI, fieldJ;
5145: PetscInt totDim, totDimAux, cStart, cEnd, numCells, c;
5146: PetscBool isMatIS, isMatISP, hasJac, hasPrec, hasDyn, hasFV = PETSC_FALSE, transform;
5147: PetscErrorCode ierr;
5150: PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);
5151: ISGetLocalSize(cellIS, &numCells);
5152: ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
5153: DMHasBasisTransform(dm, &transform);
5154: DMGetBasisTransformDM_Internal(dm, &tdm);
5155: DMGetBasisTransformVec_Internal(dm, &tv);
5156: DMGetLocalSection(dm, §ion);
5157: PetscObjectTypeCompare((PetscObject) JacP, MATIS, &isMatISP);
5158: DMGetGlobalSection(dm, &globalSection);
5159: if (isMatISP) {DMPlexGetSubdomainSection(dm, &subSection);}
5160: ISGetLocalSize(cellIS, &numCells);
5161: ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
5162: DMGetCellDS(dm, cells ? cells[cStart] : cStart, &prob);
5163: PetscDSGetNumFields(prob, &Nf);
5164: PetscDSGetTotalDimension(prob, &totDim);
5165: PetscDSHasJacobian(prob, &hasJac);
5166: PetscDSHasJacobianPreconditioner(prob, &hasPrec);
5167: /* user passed in the same matrix, avoid double contributions and
5168: only assemble the Jacobian */
5169: if (hasJac && Jac == JacP) hasPrec = PETSC_FALSE;
5170: PetscDSHasDynamicJacobian(prob, &hasDyn);
5171: hasDyn = hasDyn && (X_tShift != 0.0) ? PETSC_TRUE : PETSC_FALSE;
5172: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
5173: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
5174: if (dmAux) {
5175: DMGetEnclosureRelation(dmAux, dm, &encAux);
5176: DMConvert(dmAux, DMPLEX, &plex);
5177: DMGetLocalSection(plex, §ionAux);
5178: DMGetDS(dmAux, &probAux);
5179: PetscDSGetTotalDimension(probAux, &totDimAux);
5180: }
5181: PetscMalloc5(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,hasJac ? numCells*totDim*totDim : 0,&elemMat,hasPrec ? numCells*totDim*totDim : 0, &elemMatP,hasDyn ? numCells*totDim*totDim : 0, &elemMatD);
5182: if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
5183: DMGetCoordinateField(dm, &coordField);
5184: for (c = cStart; c < cEnd; ++c) {
5185: const PetscInt cell = cells ? cells[c] : c;
5186: const PetscInt cind = c - cStart;
5187: PetscScalar *x = NULL, *x_t = NULL;
5188: PetscInt i;
5190: DMPlexVecGetClosure(dm, section, X, cell, NULL, &x);
5191: for (i = 0; i < totDim; ++i) u[cind*totDim+i] = x[i];
5192: DMPlexVecRestoreClosure(dm, section, X, cell, NULL, &x);
5193: if (X_t) {
5194: DMPlexVecGetClosure(dm, section, X_t, cell, NULL, &x_t);
5195: for (i = 0; i < totDim; ++i) u_t[cind*totDim+i] = x_t[i];
5196: DMPlexVecRestoreClosure(dm, section, X_t, cell, NULL, &x_t);
5197: }
5198: if (dmAux) {
5199: PetscInt subcell;
5200: DMGetEnclosurePoint(dmAux, dm, encAux, cell, &subcell);
5201: DMPlexVecGetClosure(plex, sectionAux, A, subcell, NULL, &x);
5202: for (i = 0; i < totDimAux; ++i) a[cind*totDimAux+i] = x[i];
5203: DMPlexVecRestoreClosure(plex, sectionAux, A, subcell, NULL, &x);
5204: }
5205: }
5206: if (hasJac) {PetscArrayzero(elemMat, numCells*totDim*totDim);}
5207: if (hasPrec) {PetscArrayzero(elemMatP, numCells*totDim*totDim);}
5208: if (hasDyn) {PetscArrayzero(elemMatD, numCells*totDim*totDim);}
5209: for (fieldI = 0; fieldI < Nf; ++fieldI) {
5210: PetscClassId id;
5211: PetscFE fe;
5212: PetscQuadrature qGeom = NULL;
5213: PetscInt Nb;
5214: /* Conforming batches */
5215: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
5216: /* Remainder */
5217: PetscInt Nr, offset, Nq;
5218: PetscInt maxDegree;
5219: PetscFEGeom *cgeomFEM, *chunkGeom = NULL, *remGeom = NULL;
5221: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
5222: PetscObjectGetClassId((PetscObject) fe, &id);
5223: if (id == PETSCFV_CLASSID) {hasFV = PETSC_TRUE; continue;}
5224: PetscFEGetDimension(fe, &Nb);
5225: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
5226: DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);
5227: if (maxDegree <= 1) {
5228: DMFieldCreateDefaultQuadrature(coordField,cellIS,&qGeom);
5229: }
5230: if (!qGeom) {
5231: PetscFEGetQuadrature(fe,&qGeom);
5232: PetscObjectReference((PetscObject)qGeom);
5233: }
5234: PetscQuadratureGetData(qGeom, NULL, NULL, &Nq, NULL, NULL);
5235: DMSNESGetFEGeom(coordField,cellIS,qGeom,PETSC_FALSE,&cgeomFEM);
5236: blockSize = Nb;
5237: batchSize = numBlocks * blockSize;
5238: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
5239: numChunks = numCells / (numBatches*batchSize);
5240: Ne = numChunks*numBatches*batchSize;
5241: Nr = numCells % (numBatches*batchSize);
5242: offset = numCells - Nr;
5243: PetscFEGeomGetChunk(cgeomFEM,0,offset,&chunkGeom);
5244: PetscFEGeomGetChunk(cgeomFEM,offset,numCells,&remGeom);
5245: for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
5246: if (hasJac) {
5247: PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Ne, chunkGeom, u, u_t, probAux, a, t, X_tShift, elemMat);
5248: PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Nr, remGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, X_tShift, &elemMat[offset*totDim*totDim]);
5249: }
5250: if (hasPrec) {
5251: PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_PRE, fieldI, fieldJ, Ne, chunkGeom, u, u_t, probAux, a, t, X_tShift, elemMatP);
5252: PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_PRE, fieldI, fieldJ, Nr, remGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, X_tShift, &elemMatP[offset*totDim*totDim]);
5253: }
5254: if (hasDyn) {
5255: PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Ne, chunkGeom, u, u_t, probAux, a, t, X_tShift, elemMatD);
5256: PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Nr, remGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, X_tShift, &elemMatD[offset*totDim*totDim]);
5257: }
5258: }
5259: PetscFEGeomRestoreChunk(cgeomFEM,offset,numCells,&remGeom);
5260: PetscFEGeomRestoreChunk(cgeomFEM,0,offset,&chunkGeom);
5261: DMSNESRestoreFEGeom(coordField,cellIS,qGeom,PETSC_FALSE,&cgeomFEM);
5262: PetscQuadratureDestroy(&qGeom);
5263: }
5264: /* Add contribution from X_t */
5265: if (hasDyn) {for (c = 0; c < numCells*totDim*totDim; ++c) elemMat[c] += X_tShift*elemMatD[c];}
5266: if (hasFV) {
5267: PetscClassId id;
5268: PetscFV fv;
5269: PetscInt offsetI, NcI, NbI = 1, fc, f;
5271: for (fieldI = 0; fieldI < Nf; ++fieldI) {
5272: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fv);
5273: PetscDSGetFieldOffset(prob, fieldI, &offsetI);
5274: PetscObjectGetClassId((PetscObject) fv, &id);
5275: if (id != PETSCFV_CLASSID) continue;
5276: /* Put in the identity */
5277: PetscFVGetNumComponents(fv, &NcI);
5278: for (c = cStart; c < cEnd; ++c) {
5279: const PetscInt cind = c - cStart;
5280: const PetscInt eOffset = cind*totDim*totDim;
5281: for (fc = 0; fc < NcI; ++fc) {
5282: for (f = 0; f < NbI; ++f) {
5283: const PetscInt i = offsetI + f*NcI+fc;
5284: if (hasPrec) {
5285: if (hasJac) {elemMat[eOffset+i*totDim+i] = 1.0;}
5286: elemMatP[eOffset+i*totDim+i] = 1.0;
5287: } else {elemMat[eOffset+i*totDim+i] = 1.0;}
5288: }
5289: }
5290: }
5291: }
5292: /* No allocated space for FV stuff, so ignore the zero entries */
5293: MatSetOption(JacP, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);
5294: }
5295: /* Insert values into matrix */
5296: isMatIS = PETSC_FALSE;
5297: if (hasPrec && hasJac) {
5298: PetscObjectTypeCompare((PetscObject) JacP, MATIS, &isMatIS);
5299: }
5300: if (isMatIS && !subSection) {
5301: DMPlexGetSubdomainSection(dm, &subSection);
5302: }
5303: for (c = cStart; c < cEnd; ++c) {
5304: const PetscInt cell = cells ? cells[c] : c;
5305: const PetscInt cind = c - cStart;
5307: /* Transform to global basis before insertion in Jacobian */
5308: if (transform) {DMPlexBasisTransformPointTensor_Internal(dm, tdm, tv, cell, PETSC_TRUE, totDim, &elemMat[cind*totDim*totDim]);}
5309: if (hasPrec) {
5310: if (hasJac) {
5311: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, totDim, totDim, &elemMat[cind*totDim*totDim]);}
5312: if (!isMatIS) {
5313: DMPlexMatSetClosure(dm, section, globalSection, Jac, cell, &elemMat[cind*totDim*totDim], ADD_VALUES);
5314: } else {
5315: Mat lJ;
5317: MatISGetLocalMat(Jac,&lJ);
5318: DMPlexMatSetClosure(dm, section, subSection, lJ, cell, &elemMat[cind*totDim*totDim], ADD_VALUES);
5319: }
5320: }
5321: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, totDim, totDim, &elemMatP[cind*totDim*totDim]);}
5322: if (!isMatISP) {
5323: DMPlexMatSetClosure(dm, section, globalSection, JacP, cell, &elemMatP[cind*totDim*totDim], ADD_VALUES);
5324: } else {
5325: Mat lJ;
5327: MatISGetLocalMat(JacP,&lJ);
5328: DMPlexMatSetClosure(dm, section, subSection, lJ, cell, &elemMatP[cind*totDim*totDim], ADD_VALUES);
5329: }
5330: } else {
5331: if (hasJac) {
5332: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, totDim, totDim, &elemMat[cind*totDim*totDim]);}
5333: if (!isMatISP) {
5334: DMPlexMatSetClosure(dm, section, globalSection, JacP, cell, &elemMat[cind*totDim*totDim], ADD_VALUES);
5335: } else {
5336: Mat lJ;
5338: MatISGetLocalMat(JacP,&lJ);
5339: DMPlexMatSetClosure(dm, section, subSection, lJ, cell, &elemMat[cind*totDim*totDim], ADD_VALUES);
5340: }
5341: }
5342: }
5343: }
5344: ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);
5345: if (hasFV) {MatSetOption(JacP, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);}
5346: PetscFree5(u,u_t,elemMat,elemMatP,elemMatD);
5347: if (dmAux) {
5348: PetscFree(a);
5349: DMDestroy(&plex);
5350: }
5351: /* Compute boundary integrals */
5352: DMPlexComputeBdJacobian_Internal(dm, X, X_t, t, X_tShift, Jac, JacP, user);
5353: /* Assemble matrix */
5354: if (hasJac && hasPrec) {
5355: MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);
5356: MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);
5357: }
5358: MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);
5359: MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);
5360: PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);
5361: return(0);
5362: }
5364: PetscErrorCode DMPlexComputeJacobian_Hybrid_Internal(DM dm, IS cellIS, PetscReal t, PetscReal X_tShift, Vec locX, Vec locX_t, Mat Jac, Mat JacP, void *user)
5365: {
5366: DM_Plex *mesh = (DM_Plex *) dm->data;
5367: const char *name = "Hybrid Jacobian";
5368: DM dmAux = NULL;
5369: DM plex = NULL;
5370: DM plexA = NULL;
5371: DMLabel ghostLabel = NULL;
5372: PetscDS prob = NULL;
5373: PetscDS probAux = NULL;
5374: PetscSection section = NULL;
5375: DMField coordField = NULL;
5376: Vec locA;
5377: PetscScalar *u = NULL, *u_t, *a = NULL;
5378: PetscScalar *elemMat, *elemMatP;
5379: PetscSection globalSection, subSection, sectionAux;
5380: IS chunkIS;
5381: const PetscInt *cells;
5382: PetscInt *faces;
5383: PetscInt cStart, cEnd, numCells;
5384: PetscInt Nf, fieldI, fieldJ, totDim, totDimAux, numChunks, cellChunkSize, chunk;
5385: PetscInt maxDegree = PETSC_MAX_INT;
5386: PetscQuadrature affineQuad = NULL, *quads = NULL;
5387: PetscFEGeom *affineGeom = NULL, **geoms = NULL;
5388: PetscBool isMatIS = PETSC_FALSE, isMatISP = PETSC_FALSE, hasBdJac, hasBdPrec;
5389: PetscErrorCode ierr;
5392: PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);
5393: ISGetLocalSize(cellIS, &numCells);
5394: ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
5395: DMConvert(dm, DMPLEX, &plex);
5396: DMGetSection(dm, §ion);
5397: DMGetGlobalSection(dm, &globalSection);
5398: DMGetLabel(dm, "ghost", &ghostLabel);
5399: DMGetCellDS(dm, cStart, &prob);
5400: PetscDSGetNumFields(prob, &Nf);
5401: PetscDSGetTotalDimension(prob, &totDim);
5402: PetscDSHasBdJacobian(prob, &hasBdJac);
5403: PetscDSHasBdJacobianPreconditioner(prob, &hasBdPrec);
5404: PetscObjectTypeCompare((PetscObject) JacP, MATIS, &isMatISP);
5405: if (isMatISP) {DMPlexGetSubdomainSection(plex, &subSection);}
5406: if (hasBdPrec && hasBdJac) {PetscObjectTypeCompare((PetscObject) JacP, MATIS, &isMatIS);}
5407: if (isMatIS && !subSection) {DMPlexGetSubdomainSection(plex, &subSection);}
5408: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
5409: if (locA) {
5410: VecGetDM(locA, &dmAux);
5411: DMConvert(dmAux, DMPLEX, &plexA);
5412: DMGetSection(dmAux, §ionAux);
5413: DMGetCellDS(dmAux, cStart, &probAux);
5414: PetscDSGetTotalDimension(probAux, &totDimAux);
5415: }
5416: DMGetCoordinateField(dm, &coordField);
5417: DMFieldGetDegree(coordField, cellIS, NULL, &maxDegree);
5418: if (maxDegree > 1) {
5419: PetscInt f;
5420: PetscCalloc2(Nf,&quads,Nf,&geoms);
5421: for (f = 0; f < Nf; ++f) {
5422: PetscFE fe;
5424: PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
5425: if (fe) {
5426: PetscFEGetQuadrature(fe, &quads[f]);
5427: PetscObjectReference((PetscObject) quads[f]);
5428: }
5429: }
5430: }
5431: cellChunkSize = numCells;
5432: numChunks = !numCells ? 0 : PetscCeilReal(((PetscReal) numCells)/cellChunkSize);
5433: PetscCalloc1(2*cellChunkSize, &faces);
5434: ISCreateGeneral(PETSC_COMM_SELF, cellChunkSize, faces, PETSC_USE_POINTER, &chunkIS);
5435: DMPlexGetCellFields(dm, cellIS, locX, locX_t, locA, &u, &u_t, &a);
5436: DMGetWorkArray(dm, hasBdJac ? cellChunkSize*totDim*totDim : 0, MPIU_SCALAR, &elemMat);
5437: DMGetWorkArray(dm, hasBdPrec ? cellChunkSize*totDim*totDim : 0, MPIU_SCALAR, &elemMatP);
5438: for (chunk = 0; chunk < numChunks; ++chunk) {
5439: PetscInt cS = cStart+chunk*cellChunkSize, cE = PetscMin(cS+cellChunkSize, cEnd), numCells = cE - cS, c;
5441: if (hasBdJac) {PetscMemzero(elemMat, numCells*totDim*totDim * sizeof(PetscScalar));}
5442: if (hasBdPrec) {PetscMemzero(elemMatP, numCells*totDim*totDim * sizeof(PetscScalar));}
5443: /* Get faces */
5444: for (c = cS; c < cE; ++c) {
5445: const PetscInt cell = cells ? cells[c] : c;
5446: const PetscInt *cone;
5447: DMPlexGetCone(plex, cell, &cone);
5448: faces[(c-cS)*2+0] = cone[0];
5449: faces[(c-cS)*2+1] = cone[1];
5450: }
5451: ISGeneralSetIndices(chunkIS, cellChunkSize, faces, PETSC_USE_POINTER);
5452: if (maxDegree <= 1) {
5453: if (!affineQuad) {DMFieldCreateDefaultQuadrature(coordField, chunkIS, &affineQuad);}
5454: if (affineQuad) {DMSNESGetFEGeom(coordField, chunkIS, affineQuad, PETSC_TRUE, &affineGeom);}
5455: } else {
5456: PetscInt f;
5457: for (f = 0; f < Nf; ++f) {
5458: if (quads[f]) {DMSNESGetFEGeom(coordField, chunkIS, quads[f], PETSC_TRUE, &geoms[f]);}
5459: }
5460: }
5462: for (fieldI = 0; fieldI < Nf; ++fieldI) {
5463: PetscFE feI;
5464: PetscFEGeom *geom = affineGeom ? affineGeom : geoms[fieldI];
5465: PetscFEGeom *chunkGeom = NULL, *remGeom = NULL;
5466: PetscQuadrature quad = affineQuad ? affineQuad : quads[fieldI];
5467: PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset, Nq, Nb;
5469: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &feI);
5470: if (!feI) continue;
5471: PetscFEGetTileSizes(feI, NULL, &numBlocks, NULL, &numBatches);
5472: PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
5473: PetscFEGetDimension(feI, &Nb);
5474: blockSize = Nb;
5475: batchSize = numBlocks * blockSize;
5476: PetscFESetTileSizes(feI, blockSize, numBlocks, batchSize, numBatches);
5477: numChunks = numCells / (numBatches*batchSize);
5478: Ne = numChunks*numBatches*batchSize;
5479: Nr = numCells % (numBatches*batchSize);
5480: offset = numCells - Nr;
5481: PetscFEGeomGetChunk(geom,0,offset,&chunkGeom);
5482: PetscFEGeomGetChunk(geom,offset,numCells,&remGeom);
5483: for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
5484: PetscFE feJ;
5486: PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &feJ);
5487: if (!feJ) continue;
5488: if (hasBdJac) {
5489: PetscFEIntegrateHybridJacobian(prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Ne, chunkGeom, u, u_t, probAux, a, t, X_tShift, elemMat);
5490: PetscFEIntegrateHybridJacobian(prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Nr, remGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, X_tShift, &elemMat[offset*totDim*totDim]);
5491: }
5492: if (hasBdPrec) {
5493: PetscFEIntegrateHybridJacobian(prob, PETSCFE_JACOBIAN_PRE, fieldI, fieldJ, Ne, chunkGeom, u, u_t, probAux, a, t, X_tShift, elemMatP);
5494: PetscFEIntegrateHybridJacobian(prob, PETSCFE_JACOBIAN_PRE, fieldI, fieldJ, Nr, remGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, X_tShift, &elemMatP[offset*totDim*totDim]);
5495: }
5496: }
5497: PetscFEGeomRestoreChunk(geom,offset,numCells,&remGeom);
5498: PetscFEGeomRestoreChunk(geom,0,offset,&chunkGeom);
5499: }
5500: /* Insert values into matrix */
5501: for (c = cS; c < cE; ++c) {
5502: const PetscInt cell = cells ? cells[c] : c;
5503: const PetscInt cind = c - cS;
5505: if (hasBdPrec) {
5506: if (hasBdJac) {
5507: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, totDim, totDim, &elemMat[cind*totDim*totDim]);}
5508: if (!isMatIS) {
5509: DMPlexMatSetClosure(plex, section, globalSection, Jac, cell, &elemMat[cind*totDim*totDim], ADD_VALUES);
5510: } else {
5511: Mat lJ;
5513: MatISGetLocalMat(Jac,&lJ);
5514: DMPlexMatSetClosure(plex, section, subSection, lJ, cell, &elemMat[cind*totDim*totDim], ADD_VALUES);
5515: }
5516: }
5517: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, totDim, totDim, &elemMatP[cind*totDim*totDim]);}
5518: if (!isMatISP) {
5519: DMPlexMatSetClosure(plex, section, globalSection, JacP, cell, &elemMatP[cind*totDim*totDim], ADD_VALUES);
5520: } else {
5521: Mat lJ;
5523: MatISGetLocalMat(JacP,&lJ);
5524: DMPlexMatSetClosure(plex, section, subSection, lJ, cell, &elemMatP[cind*totDim*totDim], ADD_VALUES);
5525: }
5526: } else if (hasBdJac) {
5527: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, totDim, totDim, &elemMat[cind*totDim*totDim]);}
5528: if (!isMatISP) {
5529: DMPlexMatSetClosure(plex, section, globalSection, JacP, cell, &elemMat[cind*totDim*totDim], ADD_VALUES);
5530: } else {
5531: Mat lJ;
5533: MatISGetLocalMat(JacP,&lJ);
5534: DMPlexMatSetClosure(plex, section, subSection, lJ, cell, &elemMat[cind*totDim*totDim], ADD_VALUES);
5535: }
5536: }
5537: }
5538: }
5539: DMPlexRestoreCellFields(dm, cellIS, locX, locX_t, locA, &u, &u_t, &a);
5540: DMRestoreWorkArray(dm, hasBdJac ? cellChunkSize*totDim*totDim : 0, MPIU_SCALAR, &elemMat);
5541: DMRestoreWorkArray(dm, hasBdPrec ? cellChunkSize*totDim*totDim : 0, MPIU_SCALAR, &elemMatP);
5542: PetscFree(faces);
5543: ISDestroy(&chunkIS);
5544: ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);
5545: if (maxDegree <= 1) {
5546: DMSNESRestoreFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);
5547: PetscQuadratureDestroy(&affineQuad);
5548: } else {
5549: PetscInt f;
5550: for (f = 0; f < Nf; ++f) {
5551: if (geoms) {DMSNESRestoreFEGeom(coordField,cellIS,quads[f],PETSC_FALSE, &geoms[f]);}
5552: if (quads) {PetscQuadratureDestroy(&quads[f]);}
5553: }
5554: PetscFree2(quads,geoms);
5555: }
5556: if (dmAux) {DMDestroy(&plexA);}
5557: DMDestroy(&plex);
5558: /* Assemble matrix */
5559: if (hasBdJac && hasBdPrec) {
5560: MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);
5561: MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);
5562: }
5563: MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);
5564: MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);
5565: PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);
5566: return(0);
5567: }