Actual source code: plexfem.c
petsc-3.12.5 2020-03-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 Nf, PetscScalar *mode, void *ctx)
140: {
141: const PetscInt eps[3][3][3] = {{{0, 0, 0}, {0, 0, 1}, {0, -1, 0}}, {{0, 0, -1}, {0, 0, 0}, {1, 0, 0}}, {{0, 1, 0}, {-1, 0, 0}, {0, 0, 0}}};
142: PetscInt *ctxInt = (PetscInt *) ctx;
143: PetscInt dim2 = ctxInt[0];
144: PetscInt d = ctxInt[1];
145: PetscInt i, j, k = dim > 2 ? d - dim : d;
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()
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()
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, cEndInterior, 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: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
1107: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1108: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1109: for (c = cStart; c < cEnd; ++c) {
1110: PetscScalar *x = NULL;
1111: PetscReal elemDiff = 0.0;
1112: PetscInt qc = 0;
1114: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
1115: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
1117: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1118: PetscObject obj;
1119: PetscClassId id;
1120: void * const ctx = ctxs ? ctxs[field] : NULL;
1121: PetscInt Nb, Nc, q, fc;
1123: DMGetField(dm, field, NULL, &obj);
1124: PetscObjectGetClassId(obj, &id);
1125: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1126: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1127: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1128: if (debug) {
1129: char title[1024];
1130: PetscSNPrintf(title, 1023, "Solution for Field %D", field);
1131: DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
1132: }
1133: for (q = 0; q < Nq; ++q) {
1134: PetscFEGeom qgeom;
1136: qgeom.dimEmbed = fegeom.dimEmbed;
1137: qgeom.J = &fegeom.J[q*coordDim*coordDim];
1138: qgeom.invJ = &fegeom.invJ[q*coordDim*coordDim];
1139: qgeom.detJ = &fegeom.detJ[q];
1140: 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);
1141: if (transform) {
1142: gcoords = &coords[coordDim*Nq];
1143: DMPlexBasisTransformApplyReal_Internal(dm, &coords[coordDim*q], PETSC_TRUE, coordDim, &coords[coordDim*q], gcoords, dm->transformCtx);
1144: } else {
1145: gcoords = &coords[coordDim*q];
1146: }
1147: (*funcs[field])(coordDim, time, gcoords, Nc, funcVal, ctx);
1148: if (ierr) {
1149: PetscErrorCode ierr2;
1150: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
1151: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1152: ierr2 = PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2);
1153:
1154: }
1155: if (transform) {DMPlexBasisTransformApply_Internal(dm, &coords[coordDim*q], PETSC_FALSE, Nc, funcVal, funcVal, dm->transformCtx);}
1156: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], &qgeom, q, interpolant);}
1157: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
1158: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1159: for (fc = 0; fc < Nc; ++fc) {
1160: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
1161: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %D field %D,%D diff %g\n", c, field, fc, (double)(PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q]));}
1162: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q];
1163: }
1164: }
1165: fieldOffset += Nb;
1166: qc += Nc;
1167: }
1168: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1169: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %D diff %g\n", c, (double)elemDiff);}
1170: localDiff += elemDiff;
1171: }
1172: PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
1173: MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1174: *diff = PetscSqrtReal(*diff);
1175: return(0);
1176: }
1178: 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)
1179: {
1180: const PetscInt debug = ((DM_Plex*)dm->data)->printL2;
1181: DM tdm;
1182: PetscSection section;
1183: PetscQuadrature quad;
1184: Vec localX, tv;
1185: PetscScalar *funcVal, *interpolant;
1186: const PetscReal *quadWeights;
1187: PetscFEGeom fegeom;
1188: PetscReal *coords, *gcoords;
1189: PetscReal localDiff = 0.0;
1190: PetscInt dim, coordDim, qNc = 0, Nq = 0, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset;
1191: PetscBool transform;
1192: PetscErrorCode ierr;
1195: DMGetDimension(dm, &dim);
1196: DMGetCoordinateDim(dm, &coordDim);
1197: fegeom.dimEmbed = coordDim;
1198: DMGetLocalSection(dm, §ion);
1199: PetscSectionGetNumFields(section, &numFields);
1200: DMGetLocalVector(dm, &localX);
1201: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1202: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1203: DMGetBasisTransformDM_Internal(dm, &tdm);
1204: DMGetBasisTransformVec_Internal(dm, &tv);
1205: DMHasBasisTransform(dm, &transform);
1206: for (field = 0; field < numFields; ++field) {
1207: PetscFE fe;
1208: PetscInt Nc;
1210: DMGetField(dm, field, NULL, (PetscObject *) &fe);
1211: PetscFEGetQuadrature(fe, &quad);
1212: PetscFEGetNumComponents(fe, &Nc);
1213: numComponents += Nc;
1214: }
1215: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
1216: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1217: /* DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
1218: PetscMalloc6(numComponents,&funcVal,coordDim*Nq,&coords,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ,numComponents*coordDim,&interpolant,Nq,&fegeom.detJ);
1219: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1220: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1221: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1222: for (c = cStart; c < cEnd; ++c) {
1223: PetscScalar *x = NULL;
1224: PetscReal elemDiff = 0.0;
1225: PetscInt qc = 0;
1227: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
1228: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
1230: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1231: PetscFE fe;
1232: void * const ctx = ctxs ? ctxs[field] : NULL;
1233: PetscInt Nb, Nc, q, fc;
1235: DMGetField(dm, field, NULL, (PetscObject *) &fe);
1236: PetscFEGetDimension(fe, &Nb);
1237: PetscFEGetNumComponents(fe, &Nc);
1238: if (debug) {
1239: char title[1024];
1240: PetscSNPrintf(title, 1023, "Solution for Field %D", field);
1241: DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
1242: }
1243: for (q = 0; q < Nq; ++q) {
1244: PetscFEGeom qgeom;
1246: qgeom.dimEmbed = fegeom.dimEmbed;
1247: qgeom.J = &fegeom.J[q*coordDim*coordDim];
1248: qgeom.invJ = &fegeom.invJ[q*coordDim*coordDim];
1249: qgeom.detJ = &fegeom.detJ[q];
1250: 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);
1251: if (transform) {
1252: gcoords = &coords[coordDim*Nq];
1253: DMPlexBasisTransformApplyReal_Internal(dm, &coords[coordDim*q], PETSC_TRUE, coordDim, &coords[coordDim*q], gcoords, dm->transformCtx);
1254: } else {
1255: gcoords = &coords[coordDim*q];
1256: }
1257: (*funcs[field])(coordDim, time, gcoords, n, Nc, funcVal, ctx);
1258: if (ierr) {
1259: PetscErrorCode ierr2;
1260: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
1261: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1262: ierr2 = PetscFree6(funcVal,coords,fegeom.J,fegeom.invJ,interpolant,fegeom.detJ);CHKERRQ(ierr2);
1263:
1264: }
1265: if (transform) {DMPlexBasisTransformApply_Internal(dm, &coords[coordDim*q], PETSC_FALSE, Nc, funcVal, funcVal, dm->transformCtx);}
1266: PetscFEInterpolateGradient_Static(fe, &x[fieldOffset], &qgeom, q, interpolant);
1267: /* Overwrite with the dot product if the normal is given */
1268: if (n) {
1269: for (fc = 0; fc < Nc; ++fc) {
1270: PetscScalar sum = 0.0;
1271: PetscInt d;
1272: for (d = 0; d < dim; ++d) sum += interpolant[fc*dim+d]*n[d];
1273: interpolant[fc] = sum;
1274: }
1275: }
1276: for (fc = 0; fc < Nc; ++fc) {
1277: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
1278: 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]));}
1279: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q];
1280: }
1281: }
1282: fieldOffset += Nb;
1283: qc += Nc;
1284: }
1285: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1286: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %D diff %g\n", c, (double)elemDiff);}
1287: localDiff += elemDiff;
1288: }
1289: PetscFree6(funcVal,coords,fegeom.J,fegeom.invJ,interpolant,fegeom.detJ);
1290: DMRestoreLocalVector(dm, &localX);
1291: MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1292: *diff = PetscSqrtReal(*diff);
1293: return(0);
1294: }
1296: PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
1297: {
1298: const PetscInt debug = ((DM_Plex*)dm->data)->printL2;
1299: DM tdm;
1300: PetscSection section;
1301: PetscQuadrature quad;
1302: Vec localX, tv;
1303: PetscFEGeom fegeom;
1304: PetscScalar *funcVal, *interpolant;
1305: PetscReal *coords, *gcoords;
1306: PetscReal *localDiff;
1307: const PetscReal *quadPoints, *quadWeights;
1308: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
1309: PetscBool transform;
1310: PetscErrorCode ierr;
1313: DMGetDimension(dm, &dim);
1314: DMGetCoordinateDim(dm, &coordDim);
1315: DMGetLocalSection(dm, §ion);
1316: PetscSectionGetNumFields(section, &numFields);
1317: DMGetLocalVector(dm, &localX);
1318: VecSet(localX, 0.0);
1319: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1320: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1321: DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
1322: DMGetBasisTransformDM_Internal(dm, &tdm);
1323: DMGetBasisTransformVec_Internal(dm, &tv);
1324: DMHasBasisTransform(dm, &transform);
1325: for (field = 0; field < numFields; ++field) {
1326: PetscObject obj;
1327: PetscClassId id;
1328: PetscInt Nc;
1330: DMGetField(dm, field, NULL, &obj);
1331: PetscObjectGetClassId(obj, &id);
1332: if (id == PETSCFE_CLASSID) {
1333: PetscFE fe = (PetscFE) obj;
1335: PetscFEGetQuadrature(fe, &quad);
1336: PetscFEGetNumComponents(fe, &Nc);
1337: } else if (id == PETSCFV_CLASSID) {
1338: PetscFV fv = (PetscFV) obj;
1340: PetscFVGetQuadrature(fv, &quad);
1341: PetscFVGetNumComponents(fv, &Nc);
1342: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1343: numComponents += Nc;
1344: }
1345: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1346: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1347: PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,coordDim*(Nq+1),&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);
1348: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1349: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1350: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1351: for (c = cStart; c < cEnd; ++c) {
1352: PetscScalar *x = NULL;
1353: PetscInt qc = 0;
1355: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
1356: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
1358: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1359: PetscObject obj;
1360: PetscClassId id;
1361: void * const ctx = ctxs ? ctxs[field] : NULL;
1362: PetscInt Nb, Nc, q, fc;
1364: PetscReal elemDiff = 0.0;
1366: DMGetField(dm, field, NULL, &obj);
1367: PetscObjectGetClassId(obj, &id);
1368: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1369: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1370: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1371: if (debug) {
1372: char title[1024];
1373: PetscSNPrintf(title, 1023, "Solution for Field %D", field);
1374: DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
1375: }
1376: for (q = 0; q < Nq; ++q) {
1377: PetscFEGeom qgeom;
1379: qgeom.dimEmbed = fegeom.dimEmbed;
1380: qgeom.J = &fegeom.J[q*coordDim*coordDim];
1381: qgeom.invJ = &fegeom.invJ[q*coordDim*coordDim];
1382: qgeom.detJ = &fegeom.detJ[q];
1383: 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);
1384: if (transform) {
1385: gcoords = &coords[coordDim*Nq];
1386: DMPlexBasisTransformApplyReal_Internal(dm, &coords[coordDim*q], PETSC_TRUE, coordDim, &coords[coordDim*q], gcoords, dm->transformCtx);
1387: } else {
1388: gcoords = &coords[coordDim*q];
1389: }
1390: (*funcs[field])(coordDim, time, gcoords, Nc, funcVal, ctx);
1391: if (ierr) {
1392: PetscErrorCode ierr2;
1393: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
1394: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1395: ierr2 = PetscFree7(localDiff,funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2);
1396:
1397: }
1398: if (transform) {DMPlexBasisTransformApply_Internal(dm, &coords[coordDim*q], PETSC_FALSE, Nc, funcVal, funcVal, dm->transformCtx);}
1399: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], &qgeom, q, interpolant);}
1400: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
1401: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1402: for (fc = 0; fc < Nc; ++fc) {
1403: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
1404: 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]));}
1405: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q];
1406: }
1407: }
1408: fieldOffset += Nb;
1409: qc += Nc;
1410: localDiff[field] += elemDiff;
1411: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %D field %D cum diff %g\n", c, field, (double)localDiff[field]);}
1412: }
1413: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1414: }
1415: DMRestoreLocalVector(dm, &localX);
1416: MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1417: for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]);
1418: PetscFree7(localDiff,funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
1419: return(0);
1420: }
1422: /*@C
1423: 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.
1425: Collective on dm
1427: Input Parameters:
1428: + dm - The DM
1429: . time - The time
1430: . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation
1431: . ctxs - Optional array of contexts to pass to each function, or NULL.
1432: - X - The coefficient vector u_h
1434: Output Parameter:
1435: . D - A Vec which holds the difference ||u - u_h||_2 for each cell
1437: Level: developer
1439: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
1440: @*/
1441: PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D)
1442: {
1443: PetscSection section;
1444: PetscQuadrature quad;
1445: Vec localX;
1446: PetscFEGeom fegeom;
1447: PetscScalar *funcVal, *interpolant;
1448: PetscReal *coords;
1449: const PetscReal *quadPoints, *quadWeights;
1450: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
1451: PetscErrorCode ierr;
1454: VecSet(D, 0.0);
1455: DMGetDimension(dm, &dim);
1456: DMGetCoordinateDim(dm, &coordDim);
1457: DMGetLocalSection(dm, §ion);
1458: PetscSectionGetNumFields(section, &numFields);
1459: DMGetLocalVector(dm, &localX);
1460: DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
1461: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1462: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1463: for (field = 0; field < numFields; ++field) {
1464: PetscObject obj;
1465: PetscClassId id;
1466: PetscInt Nc;
1468: DMGetField(dm, field, NULL, &obj);
1469: PetscObjectGetClassId(obj, &id);
1470: if (id == PETSCFE_CLASSID) {
1471: PetscFE fe = (PetscFE) obj;
1473: PetscFEGetQuadrature(fe, &quad);
1474: PetscFEGetNumComponents(fe, &Nc);
1475: } else if (id == PETSCFV_CLASSID) {
1476: PetscFV fv = (PetscFV) obj;
1478: PetscFVGetQuadrature(fv, &quad);
1479: PetscFVGetNumComponents(fv, &Nc);
1480: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1481: numComponents += Nc;
1482: }
1483: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1484: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1485: PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);
1486: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1487: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1488: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1489: for (c = cStart; c < cEnd; ++c) {
1490: PetscScalar *x = NULL;
1491: PetscScalar elemDiff = 0.0;
1492: PetscInt qc = 0;
1494: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
1495: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
1497: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1498: PetscObject obj;
1499: PetscClassId id;
1500: void * const ctx = ctxs ? ctxs[field] : NULL;
1501: PetscInt Nb, Nc, q, fc;
1503: DMGetField(dm, field, NULL, &obj);
1504: PetscObjectGetClassId(obj, &id);
1505: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1506: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1507: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1508: if (funcs[field]) {
1509: for (q = 0; q < Nq; ++q) {
1510: PetscFEGeom qgeom;
1512: qgeom.dimEmbed = fegeom.dimEmbed;
1513: qgeom.J = &fegeom.J[q*coordDim*coordDim];
1514: qgeom.invJ = &fegeom.invJ[q*coordDim*coordDim];
1515: qgeom.detJ = &fegeom.detJ[q];
1516: 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);
1517: (*funcs[field])(coordDim, time, &coords[q*coordDim], Nc, funcVal, ctx);
1518: if (ierr) {
1519: PetscErrorCode ierr2;
1520: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
1521: ierr2 = PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2);
1522: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1523:
1524: }
1525: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], &qgeom, q, interpolant);}
1526: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
1527: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1528: for (fc = 0; fc < Nc; ++fc) {
1529: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
1530: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q];
1531: }
1532: }
1533: }
1534: fieldOffset += Nb;
1535: qc += Nc;
1536: }
1537: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1538: VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);
1539: }
1540: PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
1541: DMRestoreLocalVector(dm, &localX);
1542: VecSqrtAbs(D);
1543: return(0);
1544: }
1546: /*@C
1547: DMPlexComputeGradientClementInterpolant - This function computes the L2 projection of the cellwise gradient of a function u onto P1, and stores it in a Vec.
1549: Collective on dm
1551: Input Parameters:
1552: + dm - The DM
1553: - LocX - The coefficient vector u_h
1555: Output Parameter:
1556: . locC - A Vec which holds the Clement interpolant of the gradient
1558: Notes:
1559: Add citation to (Clement, 1975) and definition of the interpolant
1560: \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
1562: Level: developer
1564: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
1565: @*/
1566: PetscErrorCode DMPlexComputeGradientClementInterpolant(DM dm, Vec locX, Vec locC)
1567: {
1568: DM_Plex *mesh = (DM_Plex *) dm->data;
1569: PetscInt debug = mesh->printFEM;
1570: DM dmC;
1571: PetscSection section;
1572: PetscQuadrature quad;
1573: PetscScalar *interpolant, *gradsum;
1574: PetscFEGeom fegeom;
1575: PetscReal *coords;
1576: const PetscReal *quadPoints, *quadWeights;
1577: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, vStart, vEnd, v, field, fieldOffset;
1578: PetscErrorCode ierr;
1581: VecGetDM(locC, &dmC);
1582: VecSet(locC, 0.0);
1583: DMGetDimension(dm, &dim);
1584: DMGetCoordinateDim(dm, &coordDim);
1585: fegeom.dimEmbed = coordDim;
1586: DMGetLocalSection(dm, §ion);
1587: PetscSectionGetNumFields(section, &numFields);
1588: for (field = 0; field < numFields; ++field) {
1589: PetscObject obj;
1590: PetscClassId id;
1591: PetscInt Nc;
1593: DMGetField(dm, field, NULL, &obj);
1594: PetscObjectGetClassId(obj, &id);
1595: if (id == PETSCFE_CLASSID) {
1596: PetscFE fe = (PetscFE) obj;
1598: PetscFEGetQuadrature(fe, &quad);
1599: PetscFEGetNumComponents(fe, &Nc);
1600: } else if (id == PETSCFV_CLASSID) {
1601: PetscFV fv = (PetscFV) obj;
1603: PetscFVGetQuadrature(fv, &quad);
1604: PetscFVGetNumComponents(fv, &Nc);
1605: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1606: numComponents += Nc;
1607: }
1608: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1609: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1610: PetscMalloc6(coordDim*numComponents*2,&gradsum,coordDim*numComponents,&interpolant,coordDim*Nq,&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);
1611: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1612: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1613: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1614: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1615: for (v = vStart; v < vEnd; ++v) {
1616: PetscScalar volsum = 0.0;
1617: PetscInt *star = NULL;
1618: PetscInt starSize, st, d, fc;
1620: PetscArrayzero(gradsum, coordDim*numComponents);
1621: DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1622: for (st = 0; st < starSize*2; st += 2) {
1623: const PetscInt cell = star[st];
1624: PetscScalar *grad = &gradsum[coordDim*numComponents];
1625: PetscScalar *x = NULL;
1626: PetscReal vol = 0.0;
1628: if ((cell < cStart) || (cell >= cEnd)) continue;
1629: DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
1630: DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x);
1631: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1632: PetscObject obj;
1633: PetscClassId id;
1634: PetscInt Nb, Nc, q, qc = 0;
1636: PetscArrayzero(grad, coordDim*numComponents);
1637: DMGetField(dm, field, NULL, &obj);
1638: PetscObjectGetClassId(obj, &id);
1639: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1640: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1641: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1642: for (q = 0; q < Nq; ++q) {
1643: PetscFEGeom qgeom;
1645: qgeom.dimEmbed = fegeom.dimEmbed;
1646: qgeom.J = &fegeom.J[q*coordDim*coordDim];
1647: qgeom.invJ = &fegeom.invJ[q*coordDim*coordDim];
1648: qgeom.detJ = &fegeom.detJ[q];
1649: 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);
1650: if (ierr) {
1651: PetscErrorCode ierr2;
1652: ierr2 = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr2);
1653: ierr2 = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr2);
1654: ierr2 = PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2);
1655:
1656: }
1657: if (id == PETSCFE_CLASSID) {PetscFEInterpolateGradient_Static((PetscFE) obj, &x[fieldOffset], &qgeom, q, interpolant);}
1658: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1659: for (fc = 0; fc < Nc; ++fc) {
1660: const PetscReal wt = quadWeights[q*qNc+qc+fc];
1662: for (d = 0; d < coordDim; ++d) grad[fc*coordDim+d] += interpolant[fc*dim+d]*wt*fegeom.detJ[q];
1663: }
1664: vol += quadWeights[q*qNc]*fegeom.detJ[q];
1665: }
1666: fieldOffset += Nb;
1667: qc += Nc;
1668: }
1669: DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);
1670: for (fc = 0; fc < numComponents; ++fc) {
1671: for (d = 0; d < coordDim; ++d) {
1672: gradsum[fc*coordDim+d] += grad[fc*coordDim+d];
1673: }
1674: }
1675: volsum += vol;
1676: if (debug) {
1677: PetscPrintf(PETSC_COMM_SELF, "Cell %D gradient: [", cell);
1678: for (fc = 0; fc < numComponents; ++fc) {
1679: for (d = 0; d < coordDim; ++d) {
1680: if (fc || d > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
1681: PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc*coordDim+d]));
1682: }
1683: }
1684: PetscPrintf(PETSC_COMM_SELF, "]\n");
1685: }
1686: }
1687: for (fc = 0; fc < numComponents; ++fc) {
1688: for (d = 0; d < coordDim; ++d) gradsum[fc*coordDim+d] /= volsum;
1689: }
1690: DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1691: DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES);
1692: }
1693: PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
1694: return(0);
1695: }
1697: static PetscErrorCode DMPlexComputeIntegral_Internal(DM dm, Vec X, PetscInt cStart, PetscInt cEnd, PetscScalar *cintegral, void *user)
1698: {
1699: DM dmAux = NULL;
1700: PetscDS prob, probAux = NULL;
1701: PetscSection section, sectionAux;
1702: Vec locX, locA;
1703: PetscInt dim, numCells = cEnd - cStart, c, f;
1704: PetscBool useFVM = PETSC_FALSE;
1705: /* DS */
1706: PetscInt Nf, totDim, *uOff, *uOff_x, numConstants;
1707: PetscInt NfAux, totDimAux, *aOff;
1708: PetscScalar *u, *a;
1709: const PetscScalar *constants;
1710: /* Geometry */
1711: PetscFEGeom *cgeomFEM;
1712: DM dmGrad;
1713: PetscQuadrature affineQuad = NULL;
1714: Vec cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL;
1715: PetscFVCellGeom *cgeomFVM;
1716: const PetscScalar *lgrad;
1717: PetscInt maxDegree;
1718: DMField coordField;
1719: IS cellIS;
1720: PetscErrorCode ierr;
1723: DMGetDS(dm, &prob);
1724: DMGetDimension(dm, &dim);
1725: DMGetLocalSection(dm, §ion);
1726: PetscSectionGetNumFields(section, &Nf);
1727: /* Determine which discretizations we have */
1728: for (f = 0; f < Nf; ++f) {
1729: PetscObject obj;
1730: PetscClassId id;
1732: PetscDSGetDiscretization(prob, f, &obj);
1733: PetscObjectGetClassId(obj, &id);
1734: if (id == PETSCFV_CLASSID) useFVM = PETSC_TRUE;
1735: }
1736: /* Get local solution with boundary values */
1737: DMGetLocalVector(dm, &locX);
1738: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);
1739: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
1740: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
1741: /* Read DS information */
1742: PetscDSGetTotalDimension(prob, &totDim);
1743: PetscDSGetComponentOffsets(prob, &uOff);
1744: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
1745: ISCreateStride(PETSC_COMM_SELF,numCells,cStart,1,&cellIS);
1746: PetscDSGetConstants(prob, &numConstants, &constants);
1747: /* Read Auxiliary DS information */
1748: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1749: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
1750: if (dmAux) {
1751: DMGetDS(dmAux, &probAux);
1752: PetscDSGetNumFields(probAux, &NfAux);
1753: DMGetLocalSection(dmAux, §ionAux);
1754: PetscDSGetTotalDimension(probAux, &totDimAux);
1755: PetscDSGetComponentOffsets(probAux, &aOff);
1756: }
1757: /* Allocate data arrays */
1758: PetscCalloc1(numCells*totDim, &u);
1759: if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1760: /* Read out geometry */
1761: DMGetCoordinateField(dm,&coordField);
1762: DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);
1763: if (maxDegree <= 1) {
1764: DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);
1765: if (affineQuad) {
1766: DMFieldCreateFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&cgeomFEM);
1767: }
1768: }
1769: if (useFVM) {
1770: PetscFV fv = NULL;
1771: Vec grad;
1772: PetscInt fStart, fEnd;
1773: PetscBool compGrad;
1775: for (f = 0; f < Nf; ++f) {
1776: PetscObject obj;
1777: PetscClassId id;
1779: PetscDSGetDiscretization(prob, f, &obj);
1780: PetscObjectGetClassId(obj, &id);
1781: if (id == PETSCFV_CLASSID) {fv = (PetscFV) obj; break;}
1782: }
1783: PetscFVGetComputeGradients(fv, &compGrad);
1784: PetscFVSetComputeGradients(fv, PETSC_TRUE);
1785: DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);
1786: DMPlexComputeGradientFVM(dm, fv, faceGeometryFVM, cellGeometryFVM, &dmGrad);
1787: PetscFVSetComputeGradients(fv, compGrad);
1788: VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1789: /* Reconstruct and limit cell gradients */
1790: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1791: DMGetGlobalVector(dmGrad, &grad);
1792: DMPlexReconstructGradients_Internal(dm, fv, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
1793: /* Communicate gradient values */
1794: DMGetLocalVector(dmGrad, &locGrad);
1795: DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);
1796: DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);
1797: DMRestoreGlobalVector(dmGrad, &grad);
1798: /* Handle non-essential (e.g. outflow) boundary values */
1799: DMPlexInsertBoundaryValues(dm, PETSC_FALSE, locX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);
1800: VecGetArrayRead(locGrad, &lgrad);
1801: }
1802: /* Read out data from inputs */
1803: for (c = cStart; c < cEnd; ++c) {
1804: PetscScalar *x = NULL;
1805: PetscInt i;
1807: DMPlexVecGetClosure(dm, section, locX, c, NULL, &x);
1808: for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1809: DMPlexVecRestoreClosure(dm, section, locX, c, NULL, &x);
1810: if (dmAux) {
1811: DMPlexVecGetClosure(dmAux, sectionAux, locA, c, NULL, &x);
1812: for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1813: DMPlexVecRestoreClosure(dmAux, sectionAux, locA, c, NULL, &x);
1814: }
1815: }
1816: /* Do integration for each field */
1817: for (f = 0; f < Nf; ++f) {
1818: PetscObject obj;
1819: PetscClassId id;
1820: PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;
1822: PetscDSGetDiscretization(prob, f, &obj);
1823: PetscObjectGetClassId(obj, &id);
1824: if (id == PETSCFE_CLASSID) {
1825: PetscFE fe = (PetscFE) obj;
1826: PetscQuadrature q;
1827: PetscFEGeom *chunkGeom = NULL;
1828: PetscInt Nq, Nb;
1830: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1831: PetscFEGetQuadrature(fe, &q);
1832: PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);
1833: PetscFEGetDimension(fe, &Nb);
1834: blockSize = Nb*Nq;
1835: batchSize = numBlocks * blockSize;
1836: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1837: numChunks = numCells / (numBatches*batchSize);
1838: Ne = numChunks*numBatches*batchSize;
1839: Nr = numCells % (numBatches*batchSize);
1840: offset = numCells - Nr;
1841: if (!affineQuad) {
1842: DMFieldCreateFEGeom(coordField,cellIS,q,PETSC_FALSE,&cgeomFEM);
1843: }
1844: PetscFEGeomGetChunk(cgeomFEM,0,offset,&chunkGeom);
1845: PetscFEIntegrate(prob, f, Ne, chunkGeom, u, probAux, a, cintegral);
1846: PetscFEGeomGetChunk(cgeomFEM,offset,numCells,&chunkGeom);
1847: PetscFEIntegrate(prob, f, Nr, chunkGeom, &u[offset*totDim], probAux, &a[offset*totDimAux], &cintegral[offset*Nf]);
1848: PetscFEGeomRestoreChunk(cgeomFEM,offset,numCells,&chunkGeom);
1849: if (!affineQuad) {
1850: PetscFEGeomDestroy(&cgeomFEM);
1851: }
1852: } else if (id == PETSCFV_CLASSID) {
1853: PetscInt foff;
1854: PetscPointFunc obj_func;
1855: PetscScalar lint;
1857: PetscDSGetObjective(prob, f, &obj_func);
1858: PetscDSGetFieldOffset(prob, f, &foff);
1859: if (obj_func) {
1860: for (c = 0; c < numCells; ++c) {
1861: PetscScalar *u_x;
1863: DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);
1864: 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);
1865: cintegral[c*Nf+f] += PetscRealPart(lint)*cgeomFVM[c].volume;
1866: }
1867: }
1868: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f);
1869: }
1870: /* Cleanup data arrays */
1871: if (useFVM) {
1872: VecRestoreArrayRead(locGrad, &lgrad);
1873: VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1874: DMRestoreLocalVector(dmGrad, &locGrad);
1875: VecDestroy(&faceGeometryFVM);
1876: VecDestroy(&cellGeometryFVM);
1877: DMDestroy(&dmGrad);
1878: }
1879: if (dmAux) {PetscFree(a);}
1880: PetscFree(u);
1881: /* Cleanup */
1882: if (affineQuad) {
1883: PetscFEGeomDestroy(&cgeomFEM);
1884: }
1885: PetscQuadratureDestroy(&affineQuad);
1886: ISDestroy(&cellIS);
1887: DMRestoreLocalVector(dm, &locX);
1888: return(0);
1889: }
1891: /*@
1892: DMPlexComputeIntegralFEM - Form the integral over the domain from the global input X using pointwise functions specified by the user
1894: Input Parameters:
1895: + dm - The mesh
1896: . X - Global input vector
1897: - user - The user context
1899: Output Parameter:
1900: . integral - Integral for each field
1902: Level: developer
1904: .seealso: DMPlexComputeResidualFEM()
1905: @*/
1906: PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscScalar *integral, void *user)
1907: {
1908: DM_Plex *mesh = (DM_Plex *) dm->data;
1909: PetscScalar *cintegral, *lintegral;
1910: PetscInt Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell;
1917: PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1918: DMGetNumFields(dm, &Nf);
1919: DMPlexGetVTKCellHeight(dm, &cellHeight);
1920: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
1921: DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);
1922: cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight];
1923: /* TODO Introduce a loop over large chunks (right now this is a single chunk) */
1924: PetscCalloc2(Nf, &lintegral, (cEnd-cStart)*Nf, &cintegral);
1925: DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);
1926: /* Sum up values */
1927: for (cell = cStart; cell < cEnd; ++cell) {
1928: const PetscInt c = cell - cStart;
1930: if (mesh->printFEM > 1) {DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);}
1931: for (f = 0; f < Nf; ++f) lintegral[f] += cintegral[c*Nf+f];
1932: }
1933: MPIU_Allreduce(lintegral, integral, Nf, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject) dm));
1934: if (mesh->printFEM) {
1935: PetscPrintf(PetscObjectComm((PetscObject) dm), "Integral:");
1936: for (f = 0; f < Nf; ++f) {PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", (double) PetscRealPart(integral[f]));}
1937: PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");
1938: }
1939: PetscFree2(lintegral, cintegral);
1940: PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1941: return(0);
1942: }
1944: /*@
1945: DMPlexComputeCellwiseIntegralFEM - Form the vector of cellwise integrals F from the global input X using pointwise functions specified by the user
1947: Input Parameters:
1948: + dm - The mesh
1949: . X - Global input vector
1950: - user - The user context
1952: Output Parameter:
1953: . integral - Cellwise integrals for each field
1955: Level: developer
1957: .seealso: DMPlexComputeResidualFEM()
1958: @*/
1959: PetscErrorCode DMPlexComputeCellwiseIntegralFEM(DM dm, Vec X, Vec F, void *user)
1960: {
1961: DM_Plex *mesh = (DM_Plex *) dm->data;
1962: DM dmF;
1963: PetscSection sectionF;
1964: PetscScalar *cintegral, *af;
1965: PetscInt Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell;
1972: PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1973: DMGetNumFields(dm, &Nf);
1974: DMPlexGetVTKCellHeight(dm, &cellHeight);
1975: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
1976: DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);
1977: cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight];
1978: /* TODO Introduce a loop over large chunks (right now this is a single chunk) */
1979: PetscCalloc1((cEnd-cStart)*Nf, &cintegral);
1980: DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);
1981: /* Put values in F*/
1982: VecGetDM(F, &dmF);
1983: DMGetLocalSection(dmF, §ionF);
1984: VecGetArray(F, &af);
1985: for (cell = cStart; cell < cEnd; ++cell) {
1986: const PetscInt c = cell - cStart;
1987: PetscInt dof, off;
1989: if (mesh->printFEM > 1) {DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);}
1990: PetscSectionGetDof(sectionF, cell, &dof);
1991: PetscSectionGetOffset(sectionF, cell, &off);
1992: if (dof != Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cell dofs %D != %D", dof, Nf);
1993: for (f = 0; f < Nf; ++f) af[off+f] = cintegral[c*Nf+f];
1994: }
1995: VecRestoreArray(F, &af);
1996: PetscFree(cintegral);
1997: PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1998: return(0);
1999: }
2001: static PetscErrorCode DMPlexComputeBdIntegral_Internal(DM dm, Vec locX, IS pointIS,
2002: void (*func)(PetscInt, PetscInt, PetscInt,
2003: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
2004: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
2005: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
2006: PetscScalar *fintegral, void *user)
2007: {
2008: DM plex = NULL, plexA = NULL;
2009: PetscDS prob, probAux = NULL;
2010: PetscSection section, sectionAux = NULL;
2011: Vec locA = NULL;
2012: DMField coordField;
2013: PetscInt Nf, totDim, *uOff, *uOff_x;
2014: PetscInt NfAux = 0, totDimAux = 0, *aOff = NULL;
2015: PetscScalar *u, *a = NULL;
2016: const PetscScalar *constants;
2017: PetscInt numConstants, f;
2018: PetscErrorCode ierr;
2021: DMGetCoordinateField(dm, &coordField);
2022: DMConvert(dm, DMPLEX, &plex);
2023: DMGetDS(dm, &prob);
2024: DMGetLocalSection(dm, §ion);
2025: PetscSectionGetNumFields(section, &Nf);
2026: /* Determine which discretizations we have */
2027: for (f = 0; f < Nf; ++f) {
2028: PetscObject obj;
2029: PetscClassId id;
2031: PetscDSGetDiscretization(prob, f, &obj);
2032: PetscObjectGetClassId(obj, &id);
2033: if (id == PETSCFV_CLASSID) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Not supported for FVM (field %D)", f);
2034: }
2035: /* Read DS information */
2036: PetscDSGetTotalDimension(prob, &totDim);
2037: PetscDSGetComponentOffsets(prob, &uOff);
2038: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
2039: PetscDSGetConstants(prob, &numConstants, &constants);
2040: /* Read Auxiliary DS information */
2041: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
2042: if (locA) {
2043: DM dmAux;
2045: VecGetDM(locA, &dmAux);
2046: DMConvert(dmAux, DMPLEX, &plexA);
2047: DMGetDS(dmAux, &probAux);
2048: PetscDSGetNumFields(probAux, &NfAux);
2049: DMGetLocalSection(dmAux, §ionAux);
2050: PetscDSGetTotalDimension(probAux, &totDimAux);
2051: PetscDSGetComponentOffsets(probAux, &aOff);
2052: }
2053: /* Integrate over points */
2054: {
2055: PetscFEGeom *fgeom, *chunkGeom = NULL;
2056: PetscInt maxDegree;
2057: PetscQuadrature qGeom = NULL;
2058: const PetscInt *points;
2059: PetscInt numFaces, face, Nq, field;
2060: PetscInt numChunks, chunkSize, chunk, Nr, offset;
2062: ISGetLocalSize(pointIS, &numFaces);
2063: ISGetIndices(pointIS, &points);
2064: PetscCalloc2(numFaces*totDim, &u, locA ? numFaces*totDimAux : 0, &a);
2065: DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree);
2066: for (field = 0; field < Nf; ++field) {
2067: PetscFE fe;
2069: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
2070: if (maxDegree <= 1) {DMFieldCreateDefaultQuadrature(coordField, pointIS, &qGeom);}
2071: if (!qGeom) {
2072: PetscFEGetFaceQuadrature(fe, &qGeom);
2073: PetscObjectReference((PetscObject) qGeom);
2074: }
2075: PetscQuadratureGetData(qGeom, NULL, NULL, &Nq, NULL, NULL);
2076: DMPlexGetFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom);
2077: for (face = 0; face < numFaces; ++face) {
2078: const PetscInt point = points[face], *support, *cone;
2079: PetscScalar *x = NULL;
2080: PetscInt i, coneSize, faceLoc;
2082: DMPlexGetSupport(dm, point, &support);
2083: DMPlexGetConeSize(dm, support[0], &coneSize);
2084: DMPlexGetCone(dm, support[0], &cone);
2085: for (faceLoc = 0; faceLoc < coneSize; ++faceLoc) if (cone[faceLoc] == point) break;
2086: if (faceLoc == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %D in cone of support[0] %D", face, support[0]);
2087: fgeom->face[face][0] = faceLoc;
2088: DMPlexVecGetClosure(plex, section, locX, support[0], NULL, &x);
2089: for (i = 0; i < totDim; ++i) u[face*totDim+i] = x[i];
2090: DMPlexVecRestoreClosure(plex, section, locX, support[0], NULL, &x);
2091: if (locA) {
2092: PetscInt subp;
2093: DMPlexGetSubpoint(plexA, support[0], &subp);
2094: DMPlexVecGetClosure(plexA, sectionAux, locA, subp, NULL, &x);
2095: for (i = 0; i < totDimAux; ++i) a[f*totDimAux+i] = x[i];
2096: DMPlexVecRestoreClosure(plexA, sectionAux, locA, subp, NULL, &x);
2097: }
2098: }
2099: /* Get blocking */
2100: {
2101: PetscQuadrature q;
2102: PetscInt numBatches, batchSize, numBlocks, blockSize;
2103: PetscInt Nq, Nb;
2105: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
2106: PetscFEGetQuadrature(fe, &q);
2107: PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);
2108: PetscFEGetDimension(fe, &Nb);
2109: blockSize = Nb*Nq;
2110: batchSize = numBlocks * blockSize;
2111: chunkSize = numBatches*batchSize;
2112: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2113: numChunks = numFaces / chunkSize;
2114: Nr = numFaces % chunkSize;
2115: offset = numFaces - Nr;
2116: }
2117: /* Do integration for each field */
2118: for (chunk = 0; chunk < numChunks; ++chunk) {
2119: PetscFEGeomGetChunk(fgeom, chunk*chunkSize, (chunk+1)*chunkSize, &chunkGeom);
2120: PetscFEIntegrateBd(prob, field, func, chunkSize, chunkGeom, u, probAux, a, fintegral);
2121: PetscFEGeomRestoreChunk(fgeom, 0, offset, &chunkGeom);
2122: }
2123: PetscFEGeomGetChunk(fgeom, offset, numFaces, &chunkGeom);
2124: PetscFEIntegrateBd(prob, field, func, Nr, chunkGeom, &u[offset*totDim], probAux, a ? &a[offset*totDimAux] : NULL, &fintegral[offset*Nf]);
2125: PetscFEGeomRestoreChunk(fgeom, offset, numFaces, &chunkGeom);
2126: /* Cleanup data arrays */
2127: DMPlexRestoreFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom);
2128: PetscQuadratureDestroy(&qGeom);
2129: PetscFree2(u, a);
2130: ISRestoreIndices(pointIS, &points);
2131: }
2132: }
2133: if (plex) {DMDestroy(&plex);}
2134: if (plexA) {DMDestroy(&plexA);}
2135: return(0);
2136: }
2138: /*@
2139: DMPlexComputeBdIntegral - Form the integral over the specified boundary from the global input X using pointwise functions specified by the user
2141: Input Parameters:
2142: + dm - The mesh
2143: . X - Global input vector
2144: . label - The boundary DMLabel
2145: . numVals - The number of label values to use, or PETSC_DETERMINE for all values
2146: . vals - The label values to use, or PETSC_NULL for all values
2147: . func = The function to integrate along the boundary
2148: - user - The user context
2150: Output Parameter:
2151: . integral - Integral for each field
2153: Level: developer
2155: .seealso: DMPlexComputeIntegralFEM(), DMPlexComputeBdResidualFEM()
2156: @*/
2157: PetscErrorCode DMPlexComputeBdIntegral(DM dm, Vec X, DMLabel label, PetscInt numVals, const PetscInt vals[],
2158: void (*func)(PetscInt, PetscInt, PetscInt,
2159: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
2160: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
2161: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
2162: PetscScalar *integral, void *user)
2163: {
2164: Vec locX;
2165: PetscSection section;
2166: DMLabel depthLabel;
2167: IS facetIS;
2168: PetscInt dim, Nf, f, v;
2177: PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
2178: DMPlexGetDepthLabel(dm, &depthLabel);
2179: DMGetDimension(dm, &dim);
2180: DMLabelGetStratumIS(depthLabel, dim-1, &facetIS);
2181: DMGetLocalSection(dm, §ion);
2182: PetscSectionGetNumFields(section, &Nf);
2183: /* Get local solution with boundary values */
2184: DMGetLocalVector(dm, &locX);
2185: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);
2186: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
2187: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
2188: /* Loop over label values */
2189: PetscArrayzero(integral, Nf);
2190: for (v = 0; v < numVals; ++v) {
2191: IS pointIS;
2192: PetscInt numFaces, face;
2193: PetscScalar *fintegral;
2195: DMLabelGetStratumIS(label, vals[v], &pointIS);
2196: if (!pointIS) continue; /* No points with that id on this process */
2197: {
2198: IS isectIS;
2200: /* TODO: Special cases of ISIntersect where it is quick to check a priori if one is a superset of the other */
2201: ISIntersect_Caching_Internal(facetIS, pointIS, &isectIS);
2202: ISDestroy(&pointIS);
2203: pointIS = isectIS;
2204: }
2205: ISGetLocalSize(pointIS, &numFaces);
2206: PetscCalloc1(numFaces*Nf, &fintegral);
2207: DMPlexComputeBdIntegral_Internal(dm, locX, pointIS, func, fintegral, user);
2208: /* Sum point contributions into integral */
2209: for (f = 0; f < Nf; ++f) for (face = 0; face < numFaces; ++face) integral[f] += fintegral[face*Nf+f];
2210: PetscFree(fintegral);
2211: ISDestroy(&pointIS);
2212: }
2213: DMRestoreLocalVector(dm, &locX);
2214: ISDestroy(&facetIS);
2215: PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
2216: return(0);
2217: }
2219: /*@
2220: DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM.
2222: Input Parameters:
2223: + dmf - The fine mesh
2224: . dmc - The coarse mesh
2225: - user - The user context
2227: Output Parameter:
2228: . In - The interpolation matrix
2230: Level: developer
2232: .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
2233: @*/
2234: PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user)
2235: {
2236: DM_Plex *mesh = (DM_Plex *) dmc->data;
2237: const char *name = "Interpolator";
2238: PetscDS prob;
2239: PetscFE *feRef;
2240: PetscFV *fvRef;
2241: PetscSection fsection, fglobalSection;
2242: PetscSection csection, cglobalSection;
2243: PetscScalar *elemMat;
2244: PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c;
2245: PetscInt cTotDim, rTotDim = 0;
2246: PetscErrorCode ierr;
2249: PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
2250: DMGetDimension(dmf, &dim);
2251: DMGetLocalSection(dmf, &fsection);
2252: DMGetGlobalSection(dmf, &fglobalSection);
2253: DMGetLocalSection(dmc, &csection);
2254: DMGetGlobalSection(dmc, &cglobalSection);
2255: PetscSectionGetNumFields(fsection, &Nf);
2256: DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
2257: DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
2258: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
2259: DMGetDS(dmf, &prob);
2260: PetscCalloc2(Nf,&feRef,Nf,&fvRef);
2261: for (f = 0; f < Nf; ++f) {
2262: PetscObject obj;
2263: PetscClassId id;
2264: PetscInt rNb = 0, Nc = 0;
2266: PetscDSGetDiscretization(prob, f, &obj);
2267: PetscObjectGetClassId(obj, &id);
2268: if (id == PETSCFE_CLASSID) {
2269: PetscFE fe = (PetscFE) obj;
2271: PetscFERefine(fe, &feRef[f]);
2272: PetscFEGetDimension(feRef[f], &rNb);
2273: PetscFEGetNumComponents(fe, &Nc);
2274: } else if (id == PETSCFV_CLASSID) {
2275: PetscFV fv = (PetscFV) obj;
2276: PetscDualSpace Q;
2278: PetscFVRefine(fv, &fvRef[f]);
2279: PetscFVGetDualSpace(fvRef[f], &Q);
2280: PetscDualSpaceGetDimension(Q, &rNb);
2281: PetscFVGetNumComponents(fv, &Nc);
2282: }
2283: rTotDim += rNb;
2284: }
2285: PetscDSGetTotalDimension(prob, &cTotDim);
2286: PetscMalloc1(rTotDim*cTotDim,&elemMat);
2287: PetscArrayzero(elemMat, rTotDim*cTotDim);
2288: for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
2289: PetscDualSpace Qref;
2290: PetscQuadrature f;
2291: const PetscReal *qpoints, *qweights;
2292: PetscReal *points;
2293: PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d;
2295: /* Compose points from all dual basis functionals */
2296: if (feRef[fieldI]) {
2297: PetscFEGetDualSpace(feRef[fieldI], &Qref);
2298: PetscFEGetNumComponents(feRef[fieldI], &Nc);
2299: } else {
2300: PetscFVGetDualSpace(fvRef[fieldI], &Qref);
2301: PetscFVGetNumComponents(fvRef[fieldI], &Nc);
2302: }
2303: PetscDualSpaceGetDimension(Qref, &fpdim);
2304: for (i = 0; i < fpdim; ++i) {
2305: PetscDualSpaceGetFunctional(Qref, i, &f);
2306: PetscQuadratureGetData(f, NULL, NULL, &Np, NULL, NULL);
2307: npoints += Np;
2308: }
2309: PetscMalloc1(npoints*dim,&points);
2310: for (i = 0, k = 0; i < fpdim; ++i) {
2311: PetscDualSpaceGetFunctional(Qref, i, &f);
2312: PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
2313: for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
2314: }
2316: for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
2317: PetscObject obj;
2318: PetscClassId id;
2319: PetscReal *B;
2320: PetscInt NcJ = 0, cpdim = 0, j, qNc;
2322: PetscDSGetDiscretization(prob, fieldJ, &obj);
2323: PetscObjectGetClassId(obj, &id);
2324: if (id == PETSCFE_CLASSID) {
2325: PetscFE fe = (PetscFE) obj;
2327: /* Evaluate basis at points */
2328: PetscFEGetNumComponents(fe, &NcJ);
2329: PetscFEGetDimension(fe, &cpdim);
2330: /* For now, fields only interpolate themselves */
2331: if (fieldI == fieldJ) {
2332: 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);
2333: PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);
2334: for (i = 0, k = 0; i < fpdim; ++i) {
2335: PetscDualSpaceGetFunctional(Qref, i, &f);
2336: PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
2337: 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);
2338: for (p = 0; p < Np; ++p, ++k) {
2339: for (j = 0; j < cpdim; ++j) {
2340: /*
2341: cTotDim: Total columns in element interpolation matrix, sum of number of dual basis functionals in each field
2342: offsetI, offsetJ: Offsets into the larger element interpolation matrix for different fields
2343: fpdim, i, cpdim, j: Dofs for fine and coarse grids, correspond to dual space basis functionals
2344: qNC, Nc, Ncj, c: Number of components in this field
2345: Np, p: Number of quad points in the fine grid functional i
2346: k: i*Np + p, overall point number for the interpolation
2347: */
2348: for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p*qNc+c];
2349: }
2350: }
2351: }
2352: PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);
2353: }
2354: } else if (id == PETSCFV_CLASSID) {
2355: PetscFV fv = (PetscFV) obj;
2357: /* Evaluate constant function at points */
2358: PetscFVGetNumComponents(fv, &NcJ);
2359: cpdim = 1;
2360: /* For now, fields only interpolate themselves */
2361: if (fieldI == fieldJ) {
2362: 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);
2363: for (i = 0, k = 0; i < fpdim; ++i) {
2364: PetscDualSpaceGetFunctional(Qref, i, &f);
2365: PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
2366: 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);
2367: for (p = 0; p < Np; ++p, ++k) {
2368: for (j = 0; j < cpdim; ++j) {
2369: for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += 1.0*qweights[p*qNc+c];
2370: }
2371: }
2372: }
2373: }
2374: }
2375: offsetJ += cpdim;
2376: }
2377: offsetI += fpdim;
2378: PetscFree(points);
2379: }
2380: if (mesh->printFEM > 1) {DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);}
2381: /* Preallocate matrix */
2382: {
2383: Mat preallocator;
2384: PetscScalar *vals;
2385: PetscInt *cellCIndices, *cellFIndices;
2386: PetscInt locRows, locCols, cell;
2388: MatGetLocalSize(In, &locRows, &locCols);
2389: MatCreate(PetscObjectComm((PetscObject) In), &preallocator);
2390: MatSetType(preallocator, MATPREALLOCATOR);
2391: MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);
2392: MatSetUp(preallocator);
2393: PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
2394: for (cell = cStart; cell < cEnd; ++cell) {
2395: DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
2396: MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);
2397: }
2398: PetscFree3(vals,cellCIndices,cellFIndices);
2399: MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
2400: MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
2401: MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);
2402: MatDestroy(&preallocator);
2403: }
2404: /* Fill matrix */
2405: MatZeroEntries(In);
2406: for (c = cStart; c < cEnd; ++c) {
2407: DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
2408: }
2409: for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
2410: PetscFree2(feRef,fvRef);
2411: PetscFree(elemMat);
2412: MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
2413: MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
2414: if (mesh->printFEM) {
2415: PetscPrintf(PetscObjectComm((PetscObject)In), "%s:\n", name);
2416: MatChop(In, 1.0e-10);
2417: MatView(In, NULL);
2418: }
2419: PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
2420: return(0);
2421: }
2423: PetscErrorCode DMPlexComputeMassMatrixNested(DM dmc, DM dmf, Mat mass, void *user)
2424: {
2425: SETERRQ(PetscObjectComm((PetscObject) dmc), PETSC_ERR_SUP, "Laziness");
2426: }
2428: /*@
2429: DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM.
2431: Input Parameters:
2432: + dmf - The fine mesh
2433: . dmc - The coarse mesh
2434: - user - The user context
2436: Output Parameter:
2437: . In - The interpolation matrix
2439: Level: developer
2441: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
2442: @*/
2443: PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user)
2444: {
2445: DM_Plex *mesh = (DM_Plex *) dmf->data;
2446: const char *name = "Interpolator";
2447: PetscDS prob;
2448: PetscSection fsection, csection, globalFSection, globalCSection;
2449: PetscHSetIJ ht;
2450: PetscLayout rLayout;
2451: PetscInt *dnz, *onz;
2452: PetscInt locRows, rStart, rEnd;
2453: PetscReal *x, *v0, *J, *invJ, detJ;
2454: PetscReal *v0c, *Jc, *invJc, detJc;
2455: PetscScalar *elemMat;
2456: PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell;
2460: PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
2461: DMGetCoordinateDim(dmc, &dim);
2462: DMGetDS(dmc, &prob);
2463: PetscDSGetWorkspace(prob, &x, NULL, NULL, NULL, NULL);
2464: PetscDSGetNumFields(prob, &Nf);
2465: PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
2466: PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
2467: DMGetLocalSection(dmf, &fsection);
2468: DMGetGlobalSection(dmf, &globalFSection);
2469: DMGetLocalSection(dmc, &csection);
2470: DMGetGlobalSection(dmc, &globalCSection);
2471: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
2472: PetscDSGetTotalDimension(prob, &totDim);
2473: PetscMalloc1(totDim, &elemMat);
2475: MatGetLocalSize(In, &locRows, NULL);
2476: PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
2477: PetscLayoutSetLocalSize(rLayout, locRows);
2478: PetscLayoutSetBlockSize(rLayout, 1);
2479: PetscLayoutSetUp(rLayout);
2480: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
2481: PetscLayoutDestroy(&rLayout);
2482: PetscCalloc2(locRows,&dnz,locRows,&onz);
2483: PetscHSetIJCreate(&ht);
2484: for (field = 0; field < Nf; ++field) {
2485: PetscObject obj;
2486: PetscClassId id;
2487: PetscDualSpace Q = NULL;
2488: PetscQuadrature f;
2489: const PetscReal *qpoints;
2490: PetscInt Nc, Np, fpdim, i, d;
2492: PetscDSGetDiscretization(prob, field, &obj);
2493: PetscObjectGetClassId(obj, &id);
2494: if (id == PETSCFE_CLASSID) {
2495: PetscFE fe = (PetscFE) obj;
2497: PetscFEGetDualSpace(fe, &Q);
2498: PetscFEGetNumComponents(fe, &Nc);
2499: } else if (id == PETSCFV_CLASSID) {
2500: PetscFV fv = (PetscFV) obj;
2502: PetscFVGetDualSpace(fv, &Q);
2503: Nc = 1;
2504: }
2505: PetscDualSpaceGetDimension(Q, &fpdim);
2506: /* For each fine grid cell */
2507: for (cell = cStart; cell < cEnd; ++cell) {
2508: PetscInt *findices, *cindices;
2509: PetscInt numFIndices, numCIndices;
2511: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2512: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2513: if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %D != %D dual basis vecs", numFIndices, fpdim);
2514: for (i = 0; i < fpdim; ++i) {
2515: Vec pointVec;
2516: PetscScalar *pV;
2517: PetscSF coarseCellSF = NULL;
2518: const PetscSFNode *coarseCells;
2519: PetscInt numCoarseCells, q, c;
2521: /* Get points from the dual basis functional quadrature */
2522: PetscDualSpaceGetFunctional(Q, i, &f);
2523: PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
2524: VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
2525: VecSetBlockSize(pointVec, dim);
2526: VecGetArray(pointVec, &pV);
2527: for (q = 0; q < Np; ++q) {
2528: const PetscReal xi0[3] = {-1., -1., -1.};
2530: /* Transform point to real space */
2531: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2532: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2533: }
2534: VecRestoreArray(pointVec, &pV);
2535: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2536: /* OPT: Pack all quad points from fine cell */
2537: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2538: PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
2539: /* Update preallocation info */
2540: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2541: if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2542: {
2543: PetscHashIJKey key;
2544: PetscBool missing;
2546: key.i = findices[i];
2547: if (key.i >= 0) {
2548: /* Get indices for coarse elements */
2549: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2550: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2551: for (c = 0; c < numCIndices; ++c) {
2552: key.j = cindices[c];
2553: if (key.j < 0) continue;
2554: PetscHSetIJQueryAdd(ht, key, &missing);
2555: if (missing) {
2556: if ((key.j >= rStart) && (key.j < rEnd)) ++dnz[key.i-rStart];
2557: else ++onz[key.i-rStart];
2558: }
2559: }
2560: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2561: }
2562: }
2563: }
2564: PetscSFDestroy(&coarseCellSF);
2565: VecDestroy(&pointVec);
2566: }
2567: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2568: }
2569: }
2570: PetscHSetIJDestroy(&ht);
2571: MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);
2572: MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
2573: PetscFree2(dnz,onz);
2574: for (field = 0; field < Nf; ++field) {
2575: PetscObject obj;
2576: PetscClassId id;
2577: PetscDualSpace Q = NULL;
2578: PetscQuadrature f;
2579: const PetscReal *qpoints, *qweights;
2580: PetscInt Nc, qNc, Np, fpdim, i, d;
2582: PetscDSGetDiscretization(prob, field, &obj);
2583: PetscObjectGetClassId(obj, &id);
2584: if (id == PETSCFE_CLASSID) {
2585: PetscFE fe = (PetscFE) obj;
2587: PetscFEGetDualSpace(fe, &Q);
2588: PetscFEGetNumComponents(fe, &Nc);
2589: } else if (id == PETSCFV_CLASSID) {
2590: PetscFV fv = (PetscFV) obj;
2592: PetscFVGetDualSpace(fv, &Q);
2593: Nc = 1;
2594: } else SETERRQ1(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONG,"Unknown discretization type for field %D",field);
2595: PetscDualSpaceGetDimension(Q, &fpdim);
2596: /* For each fine grid cell */
2597: for (cell = cStart; cell < cEnd; ++cell) {
2598: PetscInt *findices, *cindices;
2599: PetscInt numFIndices, numCIndices;
2601: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2602: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2603: if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %D != %D dual basis vecs", numFIndices, fpdim);
2604: for (i = 0; i < fpdim; ++i) {
2605: Vec pointVec;
2606: PetscScalar *pV;
2607: PetscSF coarseCellSF = NULL;
2608: const PetscSFNode *coarseCells;
2609: PetscInt numCoarseCells, cpdim, q, c, j;
2611: /* Get points from the dual basis functional quadrature */
2612: PetscDualSpaceGetFunctional(Q, i, &f);
2613: PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);
2614: 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);
2615: VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
2616: VecSetBlockSize(pointVec, dim);
2617: VecGetArray(pointVec, &pV);
2618: for (q = 0; q < Np; ++q) {
2619: const PetscReal xi0[3] = {-1., -1., -1.};
2621: /* Transform point to real space */
2622: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2623: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2624: }
2625: VecRestoreArray(pointVec, &pV);
2626: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2627: /* OPT: Read this out from preallocation information */
2628: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2629: /* Update preallocation info */
2630: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2631: if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2632: VecGetArray(pointVec, &pV);
2633: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2634: PetscReal pVReal[3];
2635: const PetscReal xi0[3] = {-1., -1., -1.};
2637: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2638: /* Transform points from real space to coarse reference space */
2639: DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);
2640: for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]);
2641: CoordinatesRealToRef(dim, dim, xi0, v0c, invJc, pVReal, x);
2643: if (id == PETSCFE_CLASSID) {
2644: PetscFE fe = (PetscFE) obj;
2645: PetscReal *B;
2647: /* Evaluate coarse basis on contained point */
2648: PetscFEGetDimension(fe, &cpdim);
2649: PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);
2650: PetscArrayzero(elemMat, cpdim);
2651: /* Get elemMat entries by multiplying by weight */
2652: for (j = 0; j < cpdim; ++j) {
2653: for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*qweights[ccell*qNc + c];
2654: }
2655: PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);
2656: } else {
2657: cpdim = 1;
2658: for (j = 0; j < cpdim; ++j) {
2659: for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c];
2660: }
2661: }
2662: /* Update interpolator */
2663: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2664: if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2665: MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);
2666: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2667: }
2668: VecRestoreArray(pointVec, &pV);
2669: PetscSFDestroy(&coarseCellSF);
2670: VecDestroy(&pointVec);
2671: }
2672: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2673: }
2674: }
2675: PetscFree3(v0,J,invJ);
2676: PetscFree3(v0c,Jc,invJc);
2677: PetscFree(elemMat);
2678: MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
2679: MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
2680: PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
2681: return(0);
2682: }
2684: /*@
2685: DMPlexComputeMassMatrixGeneral - Form the local portion of the mass matrix M from the coarse DM to a non-nested fine DM.
2687: Input Parameters:
2688: + dmf - The fine mesh
2689: . dmc - The coarse mesh
2690: - user - The user context
2692: Output Parameter:
2693: . mass - The mass matrix
2695: Level: developer
2697: .seealso: DMPlexComputeMassMatrixNested(), DMPlexComputeInterpolatorNested(), DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
2698: @*/
2699: PetscErrorCode DMPlexComputeMassMatrixGeneral(DM dmc, DM dmf, Mat mass, void *user)
2700: {
2701: DM_Plex *mesh = (DM_Plex *) dmf->data;
2702: const char *name = "Mass Matrix";
2703: PetscDS prob;
2704: PetscSection fsection, csection, globalFSection, globalCSection;
2705: PetscHSetIJ ht;
2706: PetscLayout rLayout;
2707: PetscInt *dnz, *onz;
2708: PetscInt locRows, rStart, rEnd;
2709: PetscReal *x, *v0, *J, *invJ, detJ;
2710: PetscReal *v0c, *Jc, *invJc, detJc;
2711: PetscScalar *elemMat;
2712: PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell;
2716: DMGetCoordinateDim(dmc, &dim);
2717: DMGetDS(dmc, &prob);
2718: PetscDSGetWorkspace(prob, &x, NULL, NULL, NULL, NULL);
2719: PetscDSGetNumFields(prob, &Nf);
2720: PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
2721: PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
2722: DMGetLocalSection(dmf, &fsection);
2723: DMGetGlobalSection(dmf, &globalFSection);
2724: DMGetLocalSection(dmc, &csection);
2725: DMGetGlobalSection(dmc, &globalCSection);
2726: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
2727: PetscDSGetTotalDimension(prob, &totDim);
2728: PetscMalloc1(totDim, &elemMat);
2730: MatGetLocalSize(mass, &locRows, NULL);
2731: PetscLayoutCreate(PetscObjectComm((PetscObject) mass), &rLayout);
2732: PetscLayoutSetLocalSize(rLayout, locRows);
2733: PetscLayoutSetBlockSize(rLayout, 1);
2734: PetscLayoutSetUp(rLayout);
2735: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
2736: PetscLayoutDestroy(&rLayout);
2737: PetscCalloc2(locRows,&dnz,locRows,&onz);
2738: PetscHSetIJCreate(&ht);
2739: for (field = 0; field < Nf; ++field) {
2740: PetscObject obj;
2741: PetscClassId id;
2742: PetscQuadrature quad;
2743: const PetscReal *qpoints;
2744: PetscInt Nq, Nc, i, d;
2746: PetscDSGetDiscretization(prob, field, &obj);
2747: PetscObjectGetClassId(obj, &id);
2748: if (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, &quad);}
2749: else {PetscFVGetQuadrature((PetscFV) obj, &quad);}
2750: PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, NULL);
2751: /* For each fine grid cell */
2752: for (cell = cStart; cell < cEnd; ++cell) {
2753: Vec pointVec;
2754: PetscScalar *pV;
2755: PetscSF coarseCellSF = NULL;
2756: const PetscSFNode *coarseCells;
2757: PetscInt numCoarseCells, q, c;
2758: PetscInt *findices, *cindices;
2759: PetscInt numFIndices, numCIndices;
2761: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2762: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2763: /* Get points from the quadrature */
2764: VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);
2765: VecSetBlockSize(pointVec, dim);
2766: VecGetArray(pointVec, &pV);
2767: for (q = 0; q < Nq; ++q) {
2768: const PetscReal xi0[3] = {-1., -1., -1.};
2770: /* Transform point to real space */
2771: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2772: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2773: }
2774: VecRestoreArray(pointVec, &pV);
2775: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2776: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2777: PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
2778: /* Update preallocation info */
2779: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2780: if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2781: {
2782: PetscHashIJKey key;
2783: PetscBool missing;
2785: for (i = 0; i < numFIndices; ++i) {
2786: key.i = findices[i];
2787: if (key.i >= 0) {
2788: /* Get indices for coarse elements */
2789: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2790: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2791: for (c = 0; c < numCIndices; ++c) {
2792: key.j = cindices[c];
2793: if (key.j < 0) continue;
2794: PetscHSetIJQueryAdd(ht, key, &missing);
2795: if (missing) {
2796: if ((key.j >= rStart) && (key.j < rEnd)) ++dnz[key.i-rStart];
2797: else ++onz[key.i-rStart];
2798: }
2799: }
2800: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2801: }
2802: }
2803: }
2804: }
2805: PetscSFDestroy(&coarseCellSF);
2806: VecDestroy(&pointVec);
2807: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2808: }
2809: }
2810: PetscHSetIJDestroy(&ht);
2811: MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
2812: MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
2813: PetscFree2(dnz,onz);
2814: for (field = 0; field < Nf; ++field) {
2815: PetscObject obj;
2816: PetscClassId id;
2817: PetscQuadrature quad;
2818: PetscReal *Bfine;
2819: const PetscReal *qpoints, *qweights;
2820: PetscInt Nq, Nc, i, d;
2822: PetscDSGetDiscretization(prob, field, &obj);
2823: PetscObjectGetClassId(obj, &id);
2824: if (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, &quad);PetscFEGetDefaultTabulation((PetscFE) obj, &Bfine, NULL, NULL);}
2825: else {PetscFVGetQuadrature((PetscFV) obj, &quad);}
2826: PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, &qweights);
2827: /* For each fine grid cell */
2828: for (cell = cStart; cell < cEnd; ++cell) {
2829: Vec pointVec;
2830: PetscScalar *pV;
2831: PetscSF coarseCellSF = NULL;
2832: const PetscSFNode *coarseCells;
2833: PetscInt numCoarseCells, cpdim, q, c, j;
2834: PetscInt *findices, *cindices;
2835: PetscInt numFIndices, numCIndices;
2837: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2838: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2839: /* Get points from the quadrature */
2840: VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);
2841: VecSetBlockSize(pointVec, dim);
2842: VecGetArray(pointVec, &pV);
2843: for (q = 0; q < Nq; ++q) {
2844: const PetscReal xi0[3] = {-1., -1., -1.};
2846: /* Transform point to real space */
2847: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2848: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2849: }
2850: VecRestoreArray(pointVec, &pV);
2851: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2852: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2853: /* Update matrix */
2854: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2855: if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2856: VecGetArray(pointVec, &pV);
2857: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2858: PetscReal pVReal[3];
2859: const PetscReal xi0[3] = {-1., -1., -1.};
2862: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2863: /* Transform points from real space to coarse reference space */
2864: DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);
2865: for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]);
2866: CoordinatesRealToRef(dim, dim, xi0, v0c, invJc, pVReal, x);
2868: if (id == PETSCFE_CLASSID) {
2869: PetscFE fe = (PetscFE) obj;
2870: PetscReal *B;
2872: /* Evaluate coarse basis on contained point */
2873: PetscFEGetDimension(fe, &cpdim);
2874: PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);
2875: /* Get elemMat entries by multiplying by weight */
2876: for (i = 0; i < numFIndices; ++i) {
2877: PetscArrayzero(elemMat, cpdim);
2878: for (j = 0; j < cpdim; ++j) {
2879: for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*Bfine[(ccell*numFIndices + i)*Nc + c]*qweights[ccell*Nc + c]*detJ;
2880: }
2881: /* Update interpolator */
2882: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2883: if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2884: MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);
2885: }
2886: PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);
2887: } else {
2888: cpdim = 1;
2889: for (i = 0; i < numFIndices; ++i) {
2890: PetscArrayzero(elemMat, cpdim);
2891: for (j = 0; j < cpdim; ++j) {
2892: for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*1.0*qweights[ccell*Nc + c]*detJ;
2893: }
2894: /* Update interpolator */
2895: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2896: PetscPrintf(PETSC_COMM_SELF, "Nq: %D %D Nf: %D %D Nc: %D %D\n", ccell, Nq, i, numFIndices, j, numCIndices);
2897: if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2898: MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);
2899: }
2900: }
2901: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2902: }
2903: VecRestoreArray(pointVec, &pV);
2904: PetscSFDestroy(&coarseCellSF);
2905: VecDestroy(&pointVec);
2906: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2907: }
2908: }
2909: PetscFree3(v0,J,invJ);
2910: PetscFree3(v0c,Jc,invJc);
2911: PetscFree(elemMat);
2912: MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
2913: MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
2914: return(0);
2915: }
2917: /*@
2918: DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns
2920: Input Parameters:
2921: + dmc - The coarse mesh
2922: - dmf - The fine mesh
2923: - user - The user context
2925: Output Parameter:
2926: . sc - The mapping
2928: Level: developer
2930: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
2931: @*/
2932: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
2933: {
2934: PetscDS prob;
2935: PetscFE *feRef;
2936: PetscFV *fvRef;
2937: Vec fv, cv;
2938: IS fis, cis;
2939: PetscSection fsection, fglobalSection, csection, cglobalSection;
2940: PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
2941: PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m;
2942: PetscBool *needAvg;
2946: PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
2947: DMGetDimension(dmf, &dim);
2948: DMGetLocalSection(dmf, &fsection);
2949: DMGetGlobalSection(dmf, &fglobalSection);
2950: DMGetLocalSection(dmc, &csection);
2951: DMGetGlobalSection(dmc, &cglobalSection);
2952: PetscSectionGetNumFields(fsection, &Nf);
2953: DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
2954: DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
2955: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
2956: DMGetDS(dmc, &prob);
2957: PetscCalloc3(Nf,&feRef,Nf,&fvRef,Nf,&needAvg);
2958: for (f = 0; f < Nf; ++f) {
2959: PetscObject obj;
2960: PetscClassId id;
2961: PetscInt fNb = 0, Nc = 0;
2963: PetscDSGetDiscretization(prob, f, &obj);
2964: PetscObjectGetClassId(obj, &id);
2965: if (id == PETSCFE_CLASSID) {
2966: PetscFE fe = (PetscFE) obj;
2967: PetscSpace sp;
2968: PetscInt maxDegree;
2970: PetscFERefine(fe, &feRef[f]);
2971: PetscFEGetDimension(feRef[f], &fNb);
2972: PetscFEGetNumComponents(fe, &Nc);
2973: PetscFEGetBasisSpace(fe, &sp);
2974: PetscSpaceGetDegree(sp, NULL, &maxDegree);
2975: if (!maxDegree) needAvg[f] = PETSC_TRUE;
2976: } else if (id == PETSCFV_CLASSID) {
2977: PetscFV fv = (PetscFV) obj;
2978: PetscDualSpace Q;
2980: PetscFVRefine(fv, &fvRef[f]);
2981: PetscFVGetDualSpace(fvRef[f], &Q);
2982: PetscDualSpaceGetDimension(Q, &fNb);
2983: PetscFVGetNumComponents(fv, &Nc);
2984: needAvg[f] = PETSC_TRUE;
2985: }
2986: fTotDim += fNb;
2987: }
2988: PetscDSGetTotalDimension(prob, &cTotDim);
2989: PetscMalloc1(cTotDim,&cmap);
2990: for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
2991: PetscFE feC;
2992: PetscFV fvC;
2993: PetscDualSpace QF, QC;
2994: PetscInt order = -1, NcF, NcC, fpdim, cpdim;
2996: if (feRef[field]) {
2997: PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
2998: PetscFEGetNumComponents(feC, &NcC);
2999: PetscFEGetNumComponents(feRef[field], &NcF);
3000: PetscFEGetDualSpace(feRef[field], &QF);
3001: PetscDualSpaceGetOrder(QF, &order);
3002: PetscDualSpaceGetDimension(QF, &fpdim);
3003: PetscFEGetDualSpace(feC, &QC);
3004: PetscDualSpaceGetDimension(QC, &cpdim);
3005: } else {
3006: PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);
3007: PetscFVGetNumComponents(fvC, &NcC);
3008: PetscFVGetNumComponents(fvRef[field], &NcF);
3009: PetscFVGetDualSpace(fvRef[field], &QF);
3010: PetscDualSpaceGetDimension(QF, &fpdim);
3011: PetscFVGetDualSpace(fvC, &QC);
3012: PetscDualSpaceGetDimension(QC, &cpdim);
3013: }
3014: 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);
3015: for (c = 0; c < cpdim; ++c) {
3016: PetscQuadrature cfunc;
3017: const PetscReal *cqpoints, *cqweights;
3018: PetscInt NqcC, NpC;
3019: PetscBool found = PETSC_FALSE;
3021: PetscDualSpaceGetFunctional(QC, c, &cfunc);
3022: PetscQuadratureGetData(cfunc, NULL, &NqcC, &NpC, &cqpoints, &cqweights);
3023: 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);
3024: if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
3025: for (f = 0; f < fpdim; ++f) {
3026: PetscQuadrature ffunc;
3027: const PetscReal *fqpoints, *fqweights;
3028: PetscReal sum = 0.0;
3029: PetscInt NqcF, NpF;
3031: PetscDualSpaceGetFunctional(QF, f, &ffunc);
3032: PetscQuadratureGetData(ffunc, NULL, &NqcF, &NpF, &fqpoints, &fqweights);
3033: 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);
3034: if (NpC != NpF) continue;
3035: for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
3036: if (sum > 1.0e-9) continue;
3037: for (d = 0; d < NcC; ++d) sum += PetscAbsReal(cqweights[d]*fqweights[d]);
3038: if (sum < 1.0e-9) continue;
3039: cmap[offsetC+c] = offsetF+f;
3040: found = PETSC_TRUE;
3041: break;
3042: }
3043: if (!found) {
3044: /* TODO We really want the average here, but some asshole put VecScatter in the interface */
3045: if (fvRef[field] || (feRef[field] && order == 0)) {
3046: cmap[offsetC+c] = offsetF+0;
3047: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection");
3048: }
3049: }
3050: offsetC += cpdim;
3051: offsetF += fpdim;
3052: }
3053: for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);PetscFVDestroy(&fvRef[f]);}
3054: PetscFree3(feRef,fvRef,needAvg);
3056: DMGetGlobalVector(dmf, &fv);
3057: DMGetGlobalVector(dmc, &cv);
3058: VecGetOwnershipRange(cv, &startC, &endC);
3059: PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
3060: PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
3061: PetscMalloc1(m,&cindices);
3062: PetscMalloc1(m,&findices);
3063: for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
3064: for (c = cStart; c < cEnd; ++c) {
3065: DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
3066: for (d = 0; d < cTotDim; ++d) {
3067: if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue;
3068: 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]]);
3069: cindices[cellCIndices[d]-startC] = cellCIndices[d];
3070: findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
3071: }
3072: }
3073: PetscFree(cmap);
3074: PetscFree2(cellCIndices,cellFIndices);
3076: ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
3077: ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
3078: VecScatterCreate(cv, cis, fv, fis, sc);
3079: ISDestroy(&cis);
3080: ISDestroy(&fis);
3081: DMRestoreGlobalVector(dmf, &fv);
3082: DMRestoreGlobalVector(dmc, &cv);
3083: PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
3084: return(0);
3085: }
3087: /*@C
3088: DMPlexGetCellFields - Retrieve the field values values for a chunk of cells
3090: Input Parameters:
3091: + dm - The DM
3092: . cellIS - The cells to include
3093: . locX - A local vector with the solution fields
3094: . locX_t - A local vector with solution field time derivatives, or NULL
3095: - locA - A local vector with auxiliary fields, or NULL
3097: Output Parameters:
3098: + u - The field coefficients
3099: . u_t - The fields derivative coefficients
3100: - a - The auxiliary field coefficients
3102: Level: developer
3104: .seealso: DMPlexGetFaceFields()
3105: @*/
3106: PetscErrorCode DMPlexGetCellFields(DM dm, IS cellIS, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
3107: {
3108: DM plex, plexA = NULL;
3109: PetscSection section, sectionAux;
3110: PetscDS prob;
3111: const PetscInt *cells;
3112: PetscInt cStart, cEnd, numCells, totDim, totDimAux, c;
3113: PetscErrorCode ierr;
3123: DMPlexConvertPlex(dm, &plex, PETSC_FALSE);
3124: ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
3125: DMGetLocalSection(dm, §ion);
3126: DMGetCellDS(dm, cStart, &prob);
3127: PetscDSGetTotalDimension(prob, &totDim);
3128: if (locA) {
3129: DM dmAux;
3130: PetscDS probAux;
3132: VecGetDM(locA, &dmAux);
3133: DMPlexConvertPlex(dmAux, &plexA, PETSC_FALSE);
3134: DMGetLocalSection(dmAux, §ionAux);
3135: DMGetDS(dmAux, &probAux);
3136: PetscDSGetTotalDimension(probAux, &totDimAux);
3137: }
3138: numCells = cEnd - cStart;
3139: DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, u);
3140: if (locX_t) {DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, u_t);} else {*u_t = NULL;}
3141: if (locA) {DMGetWorkArray(dm, numCells*totDimAux, MPIU_SCALAR, a);} else {*a = NULL;}
3142: for (c = cStart; c < cEnd; ++c) {
3143: const PetscInt cell = cells ? cells[c] : c;
3144: const PetscInt cind = c - cStart;
3145: PetscScalar *x = NULL, *x_t = NULL, *ul = *u, *ul_t = *u_t, *al = *a;
3146: PetscInt i;
3148: DMPlexVecGetClosure(plex, section, locX, cell, NULL, &x);
3149: for (i = 0; i < totDim; ++i) ul[cind*totDim+i] = x[i];
3150: DMPlexVecRestoreClosure(plex, section, locX, cell, NULL, &x);
3151: if (locX_t) {
3152: DMPlexVecGetClosure(plex, section, locX_t, cell, NULL, &x_t);
3153: for (i = 0; i < totDim; ++i) ul_t[cind*totDim+i] = x_t[i];
3154: DMPlexVecRestoreClosure(plex, section, locX_t, cell, NULL, &x_t);
3155: }
3156: if (locA) {
3157: PetscInt subcell;
3158: DMPlexGetAuxiliaryPoint(plex, plexA, cell, &subcell);
3159: DMPlexVecGetClosure(plexA, sectionAux, locA, subcell, NULL, &x);
3160: for (i = 0; i < totDimAux; ++i) al[cind*totDimAux+i] = x[i];
3161: DMPlexVecRestoreClosure(plexA, sectionAux, locA, subcell, NULL, &x);
3162: }
3163: }
3164: DMDestroy(&plex);
3165: if (locA) {DMDestroy(&plexA);}
3166: ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);
3167: return(0);
3168: }
3170: /*@C
3171: DMPlexRestoreCellFields - Restore the field values values for a chunk of cells
3173: Input Parameters:
3174: + dm - The DM
3175: . cellIS - The cells to include
3176: . locX - A local vector with the solution fields
3177: . locX_t - A local vector with solution field time derivatives, or NULL
3178: - locA - A local vector with auxiliary fields, or NULL
3180: Output Parameters:
3181: + u - The field coefficients
3182: . u_t - The fields derivative coefficients
3183: - a - The auxiliary field coefficients
3185: Level: developer
3187: .seealso: DMPlexGetFaceFields()
3188: @*/
3189: PetscErrorCode DMPlexRestoreCellFields(DM dm, IS cellIS, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
3190: {
3194: DMRestoreWorkArray(dm, 0, MPIU_SCALAR, u);
3195: if (locX_t) {DMRestoreWorkArray(dm, 0, MPIU_SCALAR, u_t);}
3196: if (locA) {DMRestoreWorkArray(dm, 0, MPIU_SCALAR, a);}
3197: return(0);
3198: }
3200: /*@C
3201: DMPlexGetFaceFields - Retrieve the field values values for a chunk of faces
3203: Input Parameters:
3204: + dm - The DM
3205: . fStart - The first face to include
3206: . fEnd - The first face to exclude
3207: . locX - A local vector with the solution fields
3208: . locX_t - A local vector with solution field time derivatives, or NULL
3209: . faceGeometry - A local vector with face geometry
3210: . cellGeometry - A local vector with cell geometry
3211: - locaGrad - A local vector with field gradients, or NULL
3213: Output Parameters:
3214: + Nface - The number of faces with field values
3215: . uL - The field values at the left side of the face
3216: - uR - The field values at the right side of the face
3218: Level: developer
3220: .seealso: DMPlexGetCellFields()
3221: @*/
3222: 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)
3223: {
3224: DM dmFace, dmCell, dmGrad = NULL;
3225: PetscSection section;
3226: PetscDS prob;
3227: DMLabel ghostLabel;
3228: const PetscScalar *facegeom, *cellgeom, *x, *lgrad;
3229: PetscBool *isFE;
3230: PetscInt dim, Nf, f, Nc, numFaces = fEnd - fStart, iface, face;
3231: PetscErrorCode ierr;
3242: DMGetDimension(dm, &dim);
3243: DMGetDS(dm, &prob);
3244: DMGetLocalSection(dm, §ion);
3245: PetscDSGetNumFields(prob, &Nf);
3246: PetscDSGetTotalComponents(prob, &Nc);
3247: PetscMalloc1(Nf, &isFE);
3248: for (f = 0; f < Nf; ++f) {
3249: PetscObject obj;
3250: PetscClassId id;
3252: PetscDSGetDiscretization(prob, f, &obj);
3253: PetscObjectGetClassId(obj, &id);
3254: if (id == PETSCFE_CLASSID) {isFE[f] = PETSC_TRUE;}
3255: else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
3256: else {isFE[f] = PETSC_FALSE;}
3257: }
3258: DMGetLabel(dm, "ghost", &ghostLabel);
3259: VecGetArrayRead(locX, &x);
3260: VecGetDM(faceGeometry, &dmFace);
3261: VecGetArrayRead(faceGeometry, &facegeom);
3262: VecGetDM(cellGeometry, &dmCell);
3263: VecGetArrayRead(cellGeometry, &cellgeom);
3264: if (locGrad) {
3265: VecGetDM(locGrad, &dmGrad);
3266: VecGetArrayRead(locGrad, &lgrad);
3267: }
3268: DMGetWorkArray(dm, numFaces*Nc, MPIU_SCALAR, uL);
3269: DMGetWorkArray(dm, numFaces*Nc, MPIU_SCALAR, uR);
3270: /* Right now just eat the extra work for FE (could make a cell loop) */
3271: for (face = fStart, iface = 0; face < fEnd; ++face) {
3272: const PetscInt *cells;
3273: PetscFVFaceGeom *fg;
3274: PetscFVCellGeom *cgL, *cgR;
3275: PetscScalar *xL, *xR, *gL, *gR;
3276: PetscScalar *uLl = *uL, *uRl = *uR;
3277: PetscInt ghost, nsupp, nchild;
3279: DMLabelGetValue(ghostLabel, face, &ghost);
3280: DMPlexGetSupportSize(dm, face, &nsupp);
3281: DMPlexGetTreeChildren(dm, face, &nchild, NULL);
3282: if (ghost >= 0 || nsupp > 2 || nchild > 0) continue;
3283: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
3284: DMPlexGetSupport(dm, face, &cells);
3285: DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
3286: DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
3287: for (f = 0; f < Nf; ++f) {
3288: PetscInt off;
3290: PetscDSGetComponentOffset(prob, f, &off);
3291: if (isFE[f]) {
3292: const PetscInt *cone;
3293: PetscInt comp, coneSizeL, coneSizeR, faceLocL, faceLocR, ldof, rdof, d;
3295: xL = xR = NULL;
3296: PetscSectionGetFieldComponents(section, f, &comp);
3297: DMPlexVecGetClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
3298: DMPlexVecGetClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
3299: DMPlexGetCone(dm, cells[0], &cone);
3300: DMPlexGetConeSize(dm, cells[0], &coneSizeL);
3301: for (faceLocL = 0; faceLocL < coneSizeL; ++faceLocL) if (cone[faceLocL] == face) break;
3302: DMPlexGetCone(dm, cells[1], &cone);
3303: DMPlexGetConeSize(dm, cells[1], &coneSizeR);
3304: for (faceLocR = 0; faceLocR < coneSizeR; ++faceLocR) if (cone[faceLocR] == face) break;
3305: 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]);
3306: /* Check that FEM field has values in the right cell (sometimes its an FV ghost cell) */
3307: /* TODO: this is a hack that might not be right for nonconforming */
3308: if (faceLocL < coneSizeL) {
3309: PetscFEEvaluateFaceFields_Internal(prob, f, faceLocL, xL, &uLl[iface*Nc+off]);
3310: if (rdof == ldof && faceLocR < coneSizeR) {PetscFEEvaluateFaceFields_Internal(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);}
3311: else {for(d = 0; d < comp; ++d) uRl[iface*Nc+off+d] = uLl[iface*Nc+off+d];}
3312: }
3313: else {
3314: PetscFEEvaluateFaceFields_Internal(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);
3315: PetscSectionGetFieldComponents(section, f, &comp);
3316: for(d = 0; d < comp; ++d) uLl[iface*Nc+off+d] = uRl[iface*Nc+off+d];
3317: }
3318: DMPlexVecRestoreClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
3319: DMPlexVecRestoreClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
3320: } else {
3321: PetscFV fv;
3322: PetscInt numComp, c;
3324: PetscDSGetDiscretization(prob, f, (PetscObject *) &fv);
3325: PetscFVGetNumComponents(fv, &numComp);
3326: DMPlexPointLocalFieldRead(dm, cells[0], f, x, &xL);
3327: DMPlexPointLocalFieldRead(dm, cells[1], f, x, &xR);
3328: if (dmGrad) {
3329: PetscReal dxL[3], dxR[3];
3331: DMPlexPointLocalRead(dmGrad, cells[0], lgrad, &gL);
3332: DMPlexPointLocalRead(dmGrad, cells[1], lgrad, &gR);
3333: DMPlex_WaxpyD_Internal(dim, -1, cgL->centroid, fg->centroid, dxL);
3334: DMPlex_WaxpyD_Internal(dim, -1, cgR->centroid, fg->centroid, dxR);
3335: for (c = 0; c < numComp; ++c) {
3336: uLl[iface*Nc+off+c] = xL[c] + DMPlex_DotD_Internal(dim, &gL[c*dim], dxL);
3337: uRl[iface*Nc+off+c] = xR[c] + DMPlex_DotD_Internal(dim, &gR[c*dim], dxR);
3338: }
3339: } else {
3340: for (c = 0; c < numComp; ++c) {
3341: uLl[iface*Nc+off+c] = xL[c];
3342: uRl[iface*Nc+off+c] = xR[c];
3343: }
3344: }
3345: }
3346: }
3347: ++iface;
3348: }
3349: *Nface = iface;
3350: VecRestoreArrayRead(locX, &x);
3351: VecRestoreArrayRead(faceGeometry, &facegeom);
3352: VecRestoreArrayRead(cellGeometry, &cellgeom);
3353: if (locGrad) {
3354: VecRestoreArrayRead(locGrad, &lgrad);
3355: }
3356: PetscFree(isFE);
3357: return(0);
3358: }
3360: /*@C
3361: DMPlexRestoreFaceFields - Restore the field values values for a chunk of faces
3363: Input Parameters:
3364: + dm - The DM
3365: . fStart - The first face to include
3366: . fEnd - The first face to exclude
3367: . locX - A local vector with the solution fields
3368: . locX_t - A local vector with solution field time derivatives, or NULL
3369: . faceGeometry - A local vector with face geometry
3370: . cellGeometry - A local vector with cell geometry
3371: - locaGrad - A local vector with field gradients, or NULL
3373: Output Parameters:
3374: + Nface - The number of faces with field values
3375: . uL - The field values at the left side of the face
3376: - uR - The field values at the right side of the face
3378: Level: developer
3380: .seealso: DMPlexGetFaceFields()
3381: @*/
3382: 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)
3383: {
3387: DMRestoreWorkArray(dm, 0, MPIU_SCALAR, uL);
3388: DMRestoreWorkArray(dm, 0, MPIU_SCALAR, uR);
3389: return(0);
3390: }
3392: /*@C
3393: DMPlexGetFaceGeometry - Retrieve the geometric values for a chunk of faces
3395: Input Parameters:
3396: + dm - The DM
3397: . fStart - The first face to include
3398: . fEnd - The first face to exclude
3399: . faceGeometry - A local vector with face geometry
3400: - cellGeometry - A local vector with cell geometry
3402: Output Parameters:
3403: + Nface - The number of faces with field values
3404: . fgeom - The extract the face centroid and normal
3405: - vol - The cell volume
3407: Level: developer
3409: .seealso: DMPlexGetCellFields()
3410: @*/
3411: PetscErrorCode DMPlexGetFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscInt *Nface, PetscFVFaceGeom **fgeom, PetscReal **vol)
3412: {
3413: DM dmFace, dmCell;
3414: DMLabel ghostLabel;
3415: const PetscScalar *facegeom, *cellgeom;
3416: PetscInt dim, numFaces = fEnd - fStart, iface, face;
3417: PetscErrorCode ierr;
3425: DMGetDimension(dm, &dim);
3426: DMGetLabel(dm, "ghost", &ghostLabel);
3427: VecGetDM(faceGeometry, &dmFace);
3428: VecGetArrayRead(faceGeometry, &facegeom);
3429: VecGetDM(cellGeometry, &dmCell);
3430: VecGetArrayRead(cellGeometry, &cellgeom);
3431: PetscMalloc1(numFaces, fgeom);
3432: DMGetWorkArray(dm, numFaces*2, MPIU_SCALAR, vol);
3433: for (face = fStart, iface = 0; face < fEnd; ++face) {
3434: const PetscInt *cells;
3435: PetscFVFaceGeom *fg;
3436: PetscFVCellGeom *cgL, *cgR;
3437: PetscFVFaceGeom *fgeoml = *fgeom;
3438: PetscReal *voll = *vol;
3439: PetscInt ghost, d, nchild, nsupp;
3441: DMLabelGetValue(ghostLabel, face, &ghost);
3442: DMPlexGetSupportSize(dm, face, &nsupp);
3443: DMPlexGetTreeChildren(dm, face, &nchild, NULL);
3444: if (ghost >= 0 || nsupp > 2 || nchild > 0) continue;
3445: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
3446: DMPlexGetSupport(dm, face, &cells);
3447: DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
3448: DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
3449: for (d = 0; d < dim; ++d) {
3450: fgeoml[iface].centroid[d] = fg->centroid[d];
3451: fgeoml[iface].normal[d] = fg->normal[d];
3452: }
3453: voll[iface*2+0] = cgL->volume;
3454: voll[iface*2+1] = cgR->volume;
3455: ++iface;
3456: }
3457: *Nface = iface;
3458: VecRestoreArrayRead(faceGeometry, &facegeom);
3459: VecRestoreArrayRead(cellGeometry, &cellgeom);
3460: return(0);
3461: }
3463: /*@C
3464: DMPlexRestoreFaceGeometry - Restore the field values values for a chunk of faces
3466: Input Parameters:
3467: + dm - The DM
3468: . fStart - The first face to include
3469: . fEnd - The first face to exclude
3470: . faceGeometry - A local vector with face geometry
3471: - cellGeometry - A local vector with cell geometry
3473: Output Parameters:
3474: + Nface - The number of faces with field values
3475: . fgeom - The extract the face centroid and normal
3476: - vol - The cell volume
3478: Level: developer
3480: .seealso: DMPlexGetFaceFields()
3481: @*/
3482: PetscErrorCode DMPlexRestoreFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscInt *Nface, PetscFVFaceGeom **fgeom, PetscReal **vol)
3483: {
3487: PetscFree(*fgeom);
3488: DMRestoreWorkArray(dm, 0, MPIU_REAL, vol);
3489: return(0);
3490: }
3492: PetscErrorCode DMSNESGetFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
3493: {
3494: char composeStr[33] = {0};
3495: PetscObjectId id;
3496: PetscContainer container;
3497: PetscErrorCode ierr;
3500: PetscObjectGetId((PetscObject)quad,&id);
3501: PetscSNPrintf(composeStr, 32, "DMSNESGetFEGeom_%x\n", id);
3502: PetscObjectQuery((PetscObject) pointIS, composeStr, (PetscObject *) &container);
3503: if (container) {
3504: PetscContainerGetPointer(container, (void **) geom);
3505: } else {
3506: DMFieldCreateFEGeom(coordField, pointIS, quad, faceData, geom);
3507: PetscContainerCreate(PETSC_COMM_SELF,&container);
3508: PetscContainerSetPointer(container, (void *) *geom);
3509: PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom);
3510: PetscObjectCompose((PetscObject) pointIS, composeStr, (PetscObject) container);
3511: PetscContainerDestroy(&container);
3512: }
3513: return(0);
3514: }
3516: PetscErrorCode DMSNESRestoreFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
3517: {
3519: *geom = NULL;
3520: return(0);
3521: }
3523: PetscErrorCode DMPlexComputeResidual_Patch_Internal(DM dm, PetscSection section, IS cellIS, PetscReal t, Vec locX, Vec locX_t, Vec locF, void *user)
3524: {
3525: DM_Plex *mesh = (DM_Plex *) dm->data;
3526: const char *name = "Residual";
3527: DM dmAux = NULL;
3528: DMLabel ghostLabel = NULL;
3529: PetscDS prob = NULL;
3530: PetscDS probAux = NULL;
3531: PetscBool useFEM = PETSC_FALSE;
3532: PetscBool isImplicit = (locX_t || t == PETSC_MIN_REAL) ? PETSC_TRUE : PETSC_FALSE;
3533: DMField coordField = NULL;
3534: Vec locA;
3535: PetscScalar *u = NULL, *u_t, *a, *uL = NULL, *uR = NULL;
3536: IS chunkIS;
3537: const PetscInt *cells;
3538: PetscInt cStart, cEnd, numCells;
3539: PetscInt Nf, f, totDim, totDimAux, numChunks, cellChunkSize, chunk, fStart, fEnd;
3540: PetscInt maxDegree = PETSC_MAX_INT;
3541: PetscQuadrature affineQuad = NULL, *quads = NULL;
3542: PetscFEGeom *affineGeom = NULL, **geoms = NULL;
3543: PetscErrorCode ierr;
3546: PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);
3547: /* FEM+FVM */
3548: /* 1: Get sizes from dm and dmAux */
3549: DMGetLabel(dm, "ghost", &ghostLabel);
3550: DMGetDS(dm, &prob);
3551: PetscDSGetNumFields(prob, &Nf);
3552: PetscDSGetTotalDimension(prob, &totDim);
3553: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
3554: if (locA) {
3555: VecGetDM(locA, &dmAux);
3556: DMGetDS(dmAux, &probAux);
3557: PetscDSGetTotalDimension(probAux, &totDimAux);
3558: }
3559: /* 2: Get geometric data */
3560: for (f = 0; f < Nf; ++f) {
3561: PetscObject obj;
3562: PetscClassId id;
3563: PetscBool fimp;
3565: PetscDSGetImplicit(prob, f, &fimp);
3566: if (isImplicit != fimp) continue;
3567: PetscDSGetDiscretization(prob, f, &obj);
3568: PetscObjectGetClassId(obj, &id);
3569: if (id == PETSCFE_CLASSID) {useFEM = PETSC_TRUE;}
3570: if (id == PETSCFV_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Use of FVM with PCPATCH not yet implemented");
3571: }
3572: if (useFEM) {
3573: DMGetCoordinateField(dm, &coordField);
3574: DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);
3575: if (maxDegree <= 1) {
3576: DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);
3577: if (affineQuad) {
3578: DMSNESGetFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);
3579: }
3580: } else {
3581: PetscCalloc2(Nf,&quads,Nf,&geoms);
3582: for (f = 0; f < Nf; ++f) {
3583: PetscObject obj;
3584: PetscClassId id;
3585: PetscBool fimp;
3587: PetscDSGetImplicit(prob, f, &fimp);
3588: if (isImplicit != fimp) continue;
3589: PetscDSGetDiscretization(prob, f, &obj);
3590: PetscObjectGetClassId(obj, &id);
3591: if (id == PETSCFE_CLASSID) {
3592: PetscFE fe = (PetscFE) obj;
3594: PetscFEGetQuadrature(fe, &quads[f]);
3595: PetscObjectReference((PetscObject)quads[f]);
3596: DMSNESGetFEGeom(coordField,cellIS,quads[f],PETSC_FALSE,&geoms[f]);
3597: }
3598: }
3599: }
3600: }
3601: /* Loop over chunks */
3602: ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
3603: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
3604: if (useFEM) {ISCreate(PETSC_COMM_SELF, &chunkIS);}
3605: numCells = cEnd - cStart;
3606: numChunks = 1;
3607: cellChunkSize = numCells/numChunks;
3608: numChunks = PetscMin(1,numCells);
3609: for (chunk = 0; chunk < numChunks; ++chunk) {
3610: PetscScalar *elemVec, *fluxL = NULL, *fluxR = NULL;
3611: PetscReal *vol = NULL;
3612: PetscFVFaceGeom *fgeom = NULL;
3613: PetscInt cS = cStart+chunk*cellChunkSize, cE = PetscMin(cS+cellChunkSize, cEnd), numCells = cE - cS, c;
3614: PetscInt numFaces = 0;
3616: /* Extract field coefficients */
3617: if (useFEM) {
3618: ISGetPointSubrange(chunkIS, cS, cE, cells);
3619: DMPlexGetCellFields(dm, chunkIS, locX, locX_t, locA, &u, &u_t, &a);
3620: DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);
3621: PetscArrayzero(elemVec, numCells*totDim);
3622: }
3623: /* TODO We will interlace both our field coefficients (u, u_t, uL, uR, etc.) and our output (elemVec, fL, fR). I think this works */
3624: /* Loop over fields */
3625: for (f = 0; f < Nf; ++f) {
3626: PetscObject obj;
3627: PetscClassId id;
3628: PetscBool fimp;
3629: PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;
3631: PetscDSGetImplicit(prob, f, &fimp);
3632: if (isImplicit != fimp) continue;
3633: PetscDSGetDiscretization(prob, f, &obj);
3634: PetscObjectGetClassId(obj, &id);
3635: if (id == PETSCFE_CLASSID) {
3636: PetscFE fe = (PetscFE) obj;
3637: PetscFEGeom *geom = affineGeom ? affineGeom : geoms[f];
3638: PetscFEGeom *chunkGeom = NULL;
3639: PetscQuadrature quad = affineQuad ? affineQuad : quads[f];
3640: PetscInt Nq, Nb;
3642: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
3643: PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
3644: PetscFEGetDimension(fe, &Nb);
3645: blockSize = Nb;
3646: batchSize = numBlocks * blockSize;
3647: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
3648: numChunks = numCells / (numBatches*batchSize);
3649: Ne = numChunks*numBatches*batchSize;
3650: Nr = numCells % (numBatches*batchSize);
3651: offset = numCells - Nr;
3652: /* Integrate FE residual to get elemVec (need fields at quadrature points) */
3653: /* 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) */
3654: PetscFEGeomGetChunk(geom,0,offset,&chunkGeom);
3655: PetscFEIntegrateResidual(prob, f, Ne, chunkGeom, u, u_t, probAux, a, t, elemVec);
3656: PetscFEGeomGetChunk(geom,offset,numCells,&chunkGeom);
3657: PetscFEIntegrateResidual(prob, f, Nr, chunkGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, &elemVec[offset*totDim]);
3658: PetscFEGeomRestoreChunk(geom,offset,numCells,&chunkGeom);
3659: } else if (id == PETSCFV_CLASSID) {
3660: PetscFV fv = (PetscFV) obj;
3662: Ne = numFaces;
3663: /* Riemann solve over faces (need fields at face centroids) */
3664: /* We need to evaluate FE fields at those coordinates */
3665: PetscFVIntegrateRHSFunction(fv, prob, f, Ne, fgeom, vol, uL, uR, fluxL, fluxR);
3666: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f);
3667: }
3668: /* Loop over domain */
3669: if (useFEM) {
3670: /* Add elemVec to locX */
3671: for (c = cS; c < cE; ++c) {
3672: const PetscInt cell = cells ? cells[c] : c;
3673: const PetscInt cind = c - cStart;
3675: if (mesh->printFEM > 1) {DMPrintCellVector(cell, name, totDim, &elemVec[cind*totDim]);}
3676: if (ghostLabel) {
3677: PetscInt ghostVal;
3679: DMLabelGetValue(ghostLabel,cell,&ghostVal);
3680: if (ghostVal > 0) continue;
3681: }
3682: DMPlexVecSetClosure(dm, section, locF, cell, &elemVec[cind*totDim], ADD_ALL_VALUES);
3683: }
3684: }
3685: /* Handle time derivative */
3686: if (locX_t) {
3687: PetscScalar *x_t, *fa;
3689: VecGetArray(locF, &fa);
3690: VecGetArray(locX_t, &x_t);
3691: for (f = 0; f < Nf; ++f) {
3692: PetscFV fv;
3693: PetscObject obj;
3694: PetscClassId id;
3695: PetscInt pdim, d;
3697: PetscDSGetDiscretization(prob, f, &obj);
3698: PetscObjectGetClassId(obj, &id);
3699: if (id != PETSCFV_CLASSID) continue;
3700: fv = (PetscFV) obj;
3701: PetscFVGetNumComponents(fv, &pdim);
3702: for (c = cS; c < cE; ++c) {
3703: const PetscInt cell = cells ? cells[c] : c;
3704: PetscScalar *u_t, *r;
3706: if (ghostLabel) {
3707: PetscInt ghostVal;
3709: DMLabelGetValue(ghostLabel, cell, &ghostVal);
3710: if (ghostVal > 0) continue;
3711: }
3712: DMPlexPointLocalFieldRead(dm, cell, f, x_t, &u_t);
3713: DMPlexPointLocalFieldRef(dm, cell, f, fa, &r);
3714: for (d = 0; d < pdim; ++d) r[d] += u_t[d];
3715: }
3716: }
3717: VecRestoreArray(locX_t, &x_t);
3718: VecRestoreArray(locF, &fa);
3719: }
3720: if (useFEM) {
3721: DMPlexRestoreCellFields(dm, chunkIS, locX, locX_t, locA, &u, &u_t, &a);
3722: DMRestoreWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);
3723: }
3724: }
3725: if (useFEM) {ISDestroy(&chunkIS);}
3726: ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);
3727: /* TODO Could include boundary residual here (see DMPlexComputeResidual_Internal) */
3728: if (useFEM) {
3729: if (maxDegree <= 1) {
3730: DMSNESRestoreFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);
3731: PetscQuadratureDestroy(&affineQuad);
3732: } else {
3733: for (f = 0; f < Nf; ++f) {
3734: DMSNESRestoreFEGeom(coordField,cellIS,quads[f],PETSC_FALSE,&geoms[f]);
3735: PetscQuadratureDestroy(&quads[f]);
3736: }
3737: PetscFree2(quads,geoms);
3738: }
3739: }
3740: PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);
3741: return(0);
3742: }
3744: /*
3745: 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
3747: X - The local solution vector
3748: X_t - The local solution time derviative vector, or NULL
3749: */
3750: PetscErrorCode DMPlexComputeJacobian_Patch_Internal(DM dm, PetscSection section, PetscSection globalSection, IS cellIS,
3751: PetscReal t, PetscReal X_tShift, Vec X, Vec X_t, Mat Jac, Mat JacP, void *ctx)
3752: {
3753: DM_Plex *mesh = (DM_Plex *) dm->data;
3754: const char *name = "Jacobian", *nameP = "JacobianPre";
3755: DM dmAux = NULL;
3756: PetscDS prob, probAux = NULL;
3757: PetscSection sectionAux = NULL;
3758: Vec A;
3759: DMField coordField;
3760: PetscFEGeom *cgeomFEM;
3761: PetscQuadrature qGeom = NULL;
3762: Mat J = Jac, JP = JacP;
3763: PetscScalar *work, *u = NULL, *u_t = NULL, *a = NULL, *elemMat = NULL, *elemMatP = NULL, *elemMatD = NULL;
3764: PetscBool hasJac, hasPrec, hasDyn, assembleJac, isMatIS, isMatISP, *isFE, hasFV = PETSC_FALSE;
3765: const PetscInt *cells;
3766: PetscInt Nf, fieldI, fieldJ, maxDegree, numCells, cStart, cEnd, numChunks, chunkSize, chunk, totDim, totDimAux = 0, sz, wsz, off = 0, offCell = 0;
3767: PetscErrorCode ierr;
3770: CHKMEMQ;
3771: ISGetLocalSize(cellIS, &numCells);
3772: ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
3773: PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);
3774: DMGetDS(dm, &prob);
3775: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
3776: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
3777: if (dmAux) {
3778: DMGetLocalSection(dmAux, §ionAux);
3779: DMGetDS(dmAux, &probAux);
3780: }
3781: /* Get flags */
3782: PetscDSGetNumFields(prob, &Nf);
3783: DMGetWorkArray(dm, Nf, MPIU_BOOL, &isFE);
3784: for (fieldI = 0; fieldI < Nf; ++fieldI) {
3785: PetscObject disc;
3786: PetscClassId id;
3787: PetscDSGetDiscretization(prob, fieldI, &disc);
3788: PetscObjectGetClassId(disc, &id);
3789: if (id == PETSCFE_CLASSID) {isFE[fieldI] = PETSC_TRUE;}
3790: else if (id == PETSCFV_CLASSID) {hasFV = PETSC_TRUE; isFE[fieldI] = PETSC_FALSE;}
3791: }
3792: PetscDSHasJacobian(prob, &hasJac);
3793: PetscDSHasJacobianPreconditioner(prob, &hasPrec);
3794: PetscDSHasDynamicJacobian(prob, &hasDyn);
3795: assembleJac = hasJac && hasPrec && (Jac != JacP) ? PETSC_TRUE : PETSC_FALSE;
3796: hasDyn = hasDyn && (X_tShift != 0.0) ? PETSC_TRUE : PETSC_FALSE;
3797: PetscObjectTypeCompare((PetscObject) Jac, MATIS, &isMatIS);
3798: PetscObjectTypeCompare((PetscObject) JacP, MATIS, &isMatISP);
3799: /* Setup input data and temp arrays (should be DMGetWorkArray) */
3800: if (isMatISP || isMatISP) {DMPlexGetSubdomainSection(dm, &globalSection);}
3801: if (isMatIS) {MatISGetLocalMat(Jac, &J);}
3802: if (isMatISP) {MatISGetLocalMat(JacP, &JP);}
3803: if (hasFV) {MatSetOption(JP, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);} /* No allocated space for FV stuff, so ignore the zero entries */
3804: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
3805: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
3806: PetscDSGetTotalDimension(prob, &totDim);
3807: if (probAux) {PetscDSGetTotalDimension(probAux, &totDimAux);}
3808: CHKMEMQ;
3809: /* Compute batch sizes */
3810: if (isFE[0]) {
3811: PetscFE fe;
3812: PetscQuadrature q;
3813: PetscInt numQuadPoints, numBatches, batchSize, numBlocks, blockSize, Nb;
3815: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
3816: PetscFEGetQuadrature(fe, &q);
3817: PetscQuadratureGetData(q, NULL, NULL, &numQuadPoints, NULL, NULL);
3818: PetscFEGetDimension(fe, &Nb);
3819: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
3820: blockSize = Nb*numQuadPoints;
3821: batchSize = numBlocks * blockSize;
3822: chunkSize = numBatches * batchSize;
3823: numChunks = numCells / chunkSize + numCells % chunkSize;
3824: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
3825: } else {
3826: chunkSize = numCells;
3827: numChunks = 1;
3828: }
3829: /* Get work space */
3830: wsz = (((X?1:0) + (X_t?1:0) + (dmAux?1:0))*totDim + ((hasJac?1:0) + (hasPrec?1:0) + (hasDyn?1:0))*totDim*totDim)*chunkSize;
3831: DMGetWorkArray(dm, wsz, MPIU_SCALAR, &work);
3832: PetscArrayzero(work, wsz);
3833: off = 0;
3834: u = X ? (sz = chunkSize*totDim, off += sz, work+off-sz) : NULL;
3835: u_t = X_t ? (sz = chunkSize*totDim, off += sz, work+off-sz) : NULL;
3836: a = dmAux ? (sz = chunkSize*totDimAux, off += sz, work+off-sz) : NULL;
3837: elemMat = hasJac ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL;
3838: elemMatP = hasPrec ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL;
3839: elemMatD = hasDyn ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL;
3840: if (off != wsz) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error is workspace size %D should be %D", off, wsz);
3841: /* Setup geometry */
3842: DMGetCoordinateField(dm, &coordField);
3843: DMFieldGetDegree(coordField, cellIS, NULL, &maxDegree);
3844: if (maxDegree <= 1) {DMFieldCreateDefaultQuadrature(coordField, cellIS, &qGeom);}
3845: if (!qGeom) {
3846: PetscFE fe;
3848: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
3849: PetscFEGetQuadrature(fe, &qGeom);
3850: PetscObjectReference((PetscObject) qGeom);
3851: }
3852: DMSNESGetFEGeom(coordField, cellIS, qGeom, PETSC_FALSE, &cgeomFEM);
3853: /* Compute volume integrals */
3854: if (assembleJac) {MatZeroEntries(J);}
3855: MatZeroEntries(JP);
3856: for (chunk = 0; chunk < numChunks; ++chunk, offCell += chunkSize) {
3857: const PetscInt Ncell = PetscMin(chunkSize, numCells - offCell);
3858: PetscInt c;
3860: /* Extract values */
3861: for (c = 0; c < Ncell; ++c) {
3862: const PetscInt cell = cells ? cells[c+offCell] : c+offCell;
3863: PetscScalar *x = NULL, *x_t = NULL;
3864: PetscInt i;
3866: if (X) {
3867: DMPlexVecGetClosure(dm, section, X, cell, NULL, &x);
3868: for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
3869: DMPlexVecRestoreClosure(dm, section, X, cell, NULL, &x);
3870: }
3871: if (X_t) {
3872: DMPlexVecGetClosure(dm, section, X_t, cell, NULL, &x_t);
3873: for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
3874: DMPlexVecRestoreClosure(dm, section, X_t, cell, NULL, &x_t);
3875: }
3876: if (dmAux) {
3877: DMPlexVecGetClosure(dmAux, sectionAux, A, cell, NULL, &x);
3878: for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
3879: DMPlexVecRestoreClosure(dmAux, sectionAux, A, cell, NULL, &x);
3880: }
3881: }
3882: CHKMEMQ;
3883: for (fieldI = 0; fieldI < Nf; ++fieldI) {
3884: PetscFE fe;
3885: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
3886: for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
3887: if (hasJac) {PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMat);}
3888: if (hasPrec) {PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_PRE, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMatP);}
3889: if (hasDyn) {PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMatD);}
3890: }
3891: /* For finite volume, add the identity */
3892: if (!isFE[fieldI]) {
3893: PetscFV fv;
3894: PetscInt eOffset = 0, Nc, fc, foff;
3896: PetscDSGetFieldOffset(prob, fieldI, &foff);
3897: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fv);
3898: PetscFVGetNumComponents(fv, &Nc);
3899: for (c = 0; c < chunkSize; ++c, eOffset += totDim*totDim) {
3900: for (fc = 0; fc < Nc; ++fc) {
3901: const PetscInt i = foff + fc;
3902: if (hasJac) {elemMat [eOffset+i*totDim+i] = 1.0;}
3903: if (hasPrec) {elemMatP[eOffset+i*totDim+i] = 1.0;}
3904: }
3905: }
3906: }
3907: }
3908: CHKMEMQ;
3909: /* Add contribution from X_t */
3910: if (hasDyn) {for (c = 0; c < chunkSize*totDim*totDim; ++c) elemMat[c] += X_tShift*elemMatD[c];}
3911: /* Insert values into matrix */
3912: for (c = 0; c < Ncell; ++c) {
3913: const PetscInt cell = cells ? cells[c+offCell] : c+offCell;
3914: if (mesh->printFEM > 1) {
3915: if (hasJac) {DMPrintCellMatrix(cell, name, totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);}
3916: if (hasPrec) {DMPrintCellMatrix(cell, nameP, totDim, totDim, &elemMatP[(c-cStart)*totDim*totDim]);}
3917: }
3918: if (assembleJac) {DMPlexMatSetClosure(dm, section, globalSection, Jac, cell, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);}
3919: DMPlexMatSetClosure(dm, section, globalSection, JP, cell, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);
3920: }
3921: CHKMEMQ;
3922: }
3923: /* Cleanup */
3924: DMSNESRestoreFEGeom(coordField, cellIS, qGeom, PETSC_FALSE, &cgeomFEM);
3925: PetscQuadratureDestroy(&qGeom);
3926: if (hasFV) {MatSetOption(JacP, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);}
3927: DMRestoreWorkArray(dm, Nf, MPIU_BOOL, &isFE);
3928: 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);
3929: /* Compute boundary integrals */
3930: /* DMPlexComputeBdJacobian_Internal(dm, X, X_t, t, X_tShift, Jac, JacP, ctx); */
3931: /* Assemble matrix */
3932: if (assembleJac) {MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);}
3933: MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);
3934: PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);
3935: CHKMEMQ;
3936: return(0);
3937: }