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