Actual source code: dmfieldds.c
1: #include <petsc/private/dmfieldimpl.h>
2: #include <petsc/private/petscfeimpl.h>
3: #include <petscfe.h>
4: #include <petscdmplex.h>
5: #include <petscds.h>
7: typedef struct _n_DMField_DS {
8: PetscBool multifieldVec;
9: PetscInt height; /* Point height at which we want values and number of discretizations */
10: PetscInt fieldNum; /* Number in DS of field which we evaluate */
11: PetscObject *disc; /* Discretizations of this field at each height */
12: Vec vec; /* Field values */
13: DM dmDG; /* DM for the DG values */
14: PetscObject *discDG; /* DG Discretizations of this field at each height */
15: Vec vecDG; /* DG Field values */
16: } DMField_DS;
18: static PetscErrorCode DMFieldDestroy_DS(DMField field)
19: {
20: DMField_DS *dsfield = (DMField_DS *)field->data;
21: PetscInt i;
23: VecDestroy(&dsfield->vec);
24: for (i = 0; i < dsfield->height; i++) PetscObjectDereference(dsfield->disc[i]);
25: PetscFree(dsfield->disc);
26: VecDestroy(&dsfield->vecDG);
27: if (dsfield->discDG)
28: for (i = 0; i < dsfield->height; i++) PetscObjectDereference(dsfield->discDG[i]);
29: PetscFree(dsfield->discDG);
30: PetscFree(dsfield);
31: return 0;
32: }
34: static PetscErrorCode DMFieldView_DS(DMField field, PetscViewer viewer)
35: {
36: DMField_DS *dsfield = (DMField_DS *)field->data;
37: PetscObject disc;
38: PetscBool iascii;
40: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
41: disc = dsfield->disc[0];
42: if (iascii) {
43: PetscViewerASCIIPrintf(viewer, "PetscDS field %" PetscInt_FMT "\n", dsfield->fieldNum);
44: PetscViewerASCIIPushTab(viewer);
45: PetscObjectView(disc, viewer);
46: PetscViewerASCIIPopTab(viewer);
47: }
48: PetscViewerASCIIPushTab(viewer);
50: VecView(dsfield->vec, viewer);
51: PetscViewerASCIIPopTab(viewer);
52: return 0;
53: }
55: static PetscErrorCode DMFieldDSGetHeightDisc(DMField field, PetscInt height, PetscObject discList[], PetscObject *disc)
56: {
57: if (!discList[height]) {
58: PetscClassId id;
60: PetscObjectGetClassId(discList[0], &id);
61: if (id == PETSCFE_CLASSID) PetscFECreateHeightTrace((PetscFE)discList[0], height, (PetscFE *)&discList[height]);
62: }
63: *disc = discList[height];
64: return 0;
65: }
67: /* y[m,c] = A[m,n,c] . b[n] */
68: #define DMFieldDSdot(y, A, b, m, n, c, cast) \
69: do { \
70: PetscInt _i, _j, _k; \
71: for (_i = 0; _i < (m); _i++) { \
72: for (_k = 0; _k < (c); _k++) (y)[_i * (c) + _k] = 0.; \
73: for (_j = 0; _j < (n); _j++) { \
74: for (_k = 0; _k < (c); _k++) (y)[_i * (c) + _k] += (A)[(_i * (n) + _j) * (c) + _k] * cast((b)[_j]); \
75: } \
76: } \
77: } while (0)
79: /*
80: Since this is used for coordinates, we need to allow for the possibility that values come from multiple sections/Vecs, so that we can have DG version of the coordinates for periodicity. This reproduces DMPlexGetCellCoordinates_Internal().
81: */
82: PetscErrorCode DMFieldGetClosure_Internal(DMField field, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *values[])
83: {
84: DMField_DS *dsfield = (DMField_DS *)field->data;
85: DM fdm = dsfield->dmDG;
86: PetscSection s = NULL;
87: const PetscScalar *cvalues;
88: PetscInt pStart, pEnd;
91: *isDG = PETSC_FALSE;
92: *Nc = 0;
93: *array = NULL;
94: *values = NULL;
95: /* Check for cellwise section */
96: if (fdm) DMGetLocalSection(fdm, &s);
97: if (!s) goto cg;
98: /* Check that the cell exists in the cellwise section */
99: PetscSectionGetChart(s, &pStart, &pEnd);
100: if (cell < pStart || cell >= pEnd) goto cg;
101: /* Check for cellwise coordinates for this cell */
102: PetscSectionGetDof(s, cell, Nc);
103: if (!*Nc) goto cg;
104: /* Check for cellwise coordinates */
105: if (!dsfield->vecDG) goto cg;
106: /* Get cellwise coordinates */
107: VecGetArrayRead(dsfield->vecDG, array);
108: DMPlexPointLocalRead(fdm, cell, *array, &cvalues);
109: DMGetWorkArray(fdm, *Nc, MPIU_SCALAR, values);
110: PetscArraycpy(*values, cvalues, *Nc);
111: VecRestoreArrayRead(dsfield->vecDG, array);
112: *isDG = PETSC_TRUE;
113: return 0;
114: cg:
115: /* Use continuous values */
116: DMFieldGetDM(field, &fdm);
117: DMGetLocalSection(fdm, &s);
118: PetscSectionGetField(s, dsfield->fieldNum, &s);
119: DMPlexVecGetClosure(fdm, s, dsfield->vec, cell, Nc, values);
120: return 0;
121: }
123: PetscErrorCode DMFieldRestoreClosure_Internal(DMField field, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *values[])
124: {
125: DMField_DS *dsfield = (DMField_DS *)field->data;
126: DM fdm;
127: PetscSection s;
130: if (*isDG) {
131: DMRestoreWorkArray(dsfield->dmDG, *Nc, MPIU_SCALAR, values);
132: } else {
133: DMFieldGetDM(field, &fdm);
134: DMGetLocalSection(fdm, &s);
135: DMPlexVecRestoreClosure(fdm, s, dsfield->vec, cell, Nc, (PetscScalar **)values);
136: }
137: return 0;
138: }
140: /* TODO: Reorganize interface so that I can reuse a tabulation rather than mallocing each time */
141: static PetscErrorCode DMFieldEvaluateFE_DS(DMField field, IS pointIS, PetscQuadrature quad, PetscDataType type, void *B, void *D, void *H)
142: {
143: DMField_DS *dsfield = (DMField_DS *)field->data;
144: DM dm;
145: PetscObject disc;
146: PetscClassId classid;
147: PetscInt nq, nc, dim, meshDim, numCells;
148: PetscSection section;
149: const PetscReal *qpoints;
150: PetscBool isStride;
151: const PetscInt *points = NULL;
152: PetscInt sfirst = -1, stride = -1;
155: dm = field->dm;
156: nc = field->numComponents;
157: PetscQuadratureGetData(quad, &dim, NULL, &nq, &qpoints, NULL);
158: DMFieldDSGetHeightDisc(field, dsfield->height - 1 - dim, dsfield->disc, &disc);
159: DMGetDimension(dm, &meshDim);
160: DMGetLocalSection(dm, §ion);
161: PetscSectionGetField(section, dsfield->fieldNum, §ion);
162: PetscObjectGetClassId(disc, &classid);
163: /* TODO: batch */
164: PetscObjectTypeCompare((PetscObject)pointIS, ISSTRIDE, &isStride);
165: ISGetLocalSize(pointIS, &numCells);
166: if (isStride) ISStrideGetInfo(pointIS, &sfirst, &stride);
167: else ISGetIndices(pointIS, &points);
168: if (classid == PETSCFE_CLASSID) {
169: PetscFE fe = (PetscFE)disc;
170: PetscInt feDim, i;
171: PetscInt K = H ? 2 : (D ? 1 : (B ? 0 : -1));
172: PetscTabulation T;
174: PetscFEGetDimension(fe, &feDim);
175: PetscFECreateTabulation(fe, 1, nq, qpoints, K, &T);
176: for (i = 0; i < numCells; i++) {
177: PetscInt c = isStride ? (sfirst + i * stride) : points[i];
178: PetscInt closureSize;
179: const PetscScalar *array;
180: PetscScalar *elem = NULL;
181: PetscBool isDG;
183: DMFieldGetClosure_Internal(field, c, &isDG, &closureSize, &array, &elem);
184: if (B) {
185: /* field[c] = T[q,b,c] . coef[b], so v[c] = T[q,b,c] . coords[b] */
186: if (type == PETSC_SCALAR) {
187: PetscScalar *cB = &((PetscScalar *)B)[nc * nq * i];
189: DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, (PetscScalar));
190: } else {
191: PetscReal *cB = &((PetscReal *)B)[nc * nq * i];
193: DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, PetscRealPart);
194: }
195: }
196: if (D) {
197: if (type == PETSC_SCALAR) {
198: PetscScalar *cD = &((PetscScalar *)D)[nc * nq * dim * i];
200: DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), (PetscScalar));
201: } else {
202: PetscReal *cD = &((PetscReal *)D)[nc * nq * dim * i];
204: DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), PetscRealPart);
205: }
206: }
207: if (H) {
208: if (type == PETSC_SCALAR) {
209: PetscScalar *cH = &((PetscScalar *)H)[nc * nq * dim * dim * i];
211: DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), (PetscScalar));
212: } else {
213: PetscReal *cH = &((PetscReal *)H)[nc * nq * dim * dim * i];
215: DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), PetscRealPart);
216: }
217: }
218: DMFieldRestoreClosure_Internal(field, c, &isDG, &closureSize, &array, &elem);
219: }
220: PetscTabulationDestroy(&T);
221: } else SETERRQ(PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "Not implemented");
222: if (!isStride) ISRestoreIndices(pointIS, &points);
223: return 0;
224: }
226: static PetscErrorCode DMFieldEvaluate_DS(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
227: {
228: DMField_DS *dsfield = (DMField_DS *)field->data;
229: PetscSF cellSF = NULL;
230: const PetscSFNode *cells;
231: PetscInt c, nFound, numCells, feDim, nc;
232: const PetscInt *cellDegrees;
233: const PetscScalar *pointsArray;
234: PetscScalar *cellPoints;
235: PetscInt gatherSize, gatherMax;
236: PetscInt dim, dimR, offset;
237: MPI_Datatype pointType;
238: PetscObject cellDisc;
239: PetscFE cellFE;
240: PetscClassId discID;
241: PetscReal *coordsReal, *coordsRef;
242: PetscSection section;
243: PetscScalar *cellBs = NULL, *cellDs = NULL, *cellHs = NULL;
244: PetscReal *cellBr = NULL, *cellDr = NULL, *cellHr = NULL;
245: PetscReal *v, *J, *invJ, *detJ;
247: nc = field->numComponents;
248: DMGetLocalSection(field->dm, §ion);
249: DMFieldDSGetHeightDisc(field, 0, dsfield->disc, &cellDisc);
250: PetscObjectGetClassId(cellDisc, &discID);
252: cellFE = (PetscFE)cellDisc;
253: PetscFEGetDimension(cellFE, &feDim);
254: DMGetCoordinateDim(field->dm, &dim);
255: DMGetDimension(field->dm, &dimR);
256: DMLocatePoints(field->dm, points, DM_POINTLOCATION_NONE, &cellSF);
257: PetscSFGetGraph(cellSF, &numCells, &nFound, NULL, &cells);
259: PetscSFComputeDegreeBegin(cellSF, &cellDegrees);
260: PetscSFComputeDegreeEnd(cellSF, &cellDegrees);
261: for (c = 0, gatherSize = 0, gatherMax = 0; c < numCells; c++) {
262: gatherMax = PetscMax(gatherMax, cellDegrees[c]);
263: gatherSize += cellDegrees[c];
264: }
265: PetscMalloc3(gatherSize * dim, &cellPoints, gatherMax * dim, &coordsReal, gatherMax * dimR, &coordsRef);
266: PetscMalloc4(gatherMax * dimR, &v, gatherMax * dimR * dimR, &J, gatherMax * dimR * dimR, &invJ, gatherMax, &detJ);
267: if (datatype == PETSC_SCALAR) PetscMalloc3(B ? nc * gatherSize : 0, &cellBs, D ? nc * dim * gatherSize : 0, &cellDs, H ? nc * dim * dim * gatherSize : 0, &cellHs);
268: else PetscMalloc3(B ? nc * gatherSize : 0, &cellBr, D ? nc * dim * gatherSize : 0, &cellDr, H ? nc * dim * dim * gatherSize : 0, &cellHr);
270: MPI_Type_contiguous(dim, MPIU_SCALAR, &pointType);
271: MPI_Type_commit(&pointType);
272: VecGetArrayRead(points, &pointsArray);
273: PetscSFGatherBegin(cellSF, pointType, pointsArray, cellPoints);
274: PetscSFGatherEnd(cellSF, pointType, pointsArray, cellPoints);
275: VecRestoreArrayRead(points, &pointsArray);
276: for (c = 0, offset = 0; c < numCells; c++) {
277: PetscInt nq = cellDegrees[c], p;
279: if (nq) {
280: PetscInt K = H ? 2 : (D ? 1 : (B ? 0 : -1));
281: PetscTabulation T;
282: PetscQuadrature quad;
283: const PetscScalar *array;
284: PetscScalar *elem = NULL;
285: PetscReal *quadPoints;
286: PetscBool isDG;
287: PetscInt closureSize, d, e, f, g;
289: for (p = 0; p < dim * nq; p++) coordsReal[p] = PetscRealPart(cellPoints[dim * offset + p]);
290: DMPlexCoordinatesToReference(field->dm, c, nq, coordsReal, coordsRef);
291: PetscFECreateTabulation(cellFE, 1, nq, coordsRef, K, &T);
292: PetscQuadratureCreate(PETSC_COMM_SELF, &quad);
293: PetscMalloc1(dimR * nq, &quadPoints);
294: for (p = 0; p < dimR * nq; p++) quadPoints[p] = coordsRef[p];
295: PetscQuadratureSetData(quad, dimR, 0, nq, quadPoints, NULL);
296: DMPlexComputeCellGeometryFEM(field->dm, c, quad, v, J, invJ, detJ);
297: PetscQuadratureDestroy(&quad);
298: DMFieldGetClosure_Internal(field, c, &isDG, &closureSize, &array, &elem);
299: if (B) {
300: if (datatype == PETSC_SCALAR) {
301: PetscScalar *cB = &cellBs[nc * offset];
303: DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, (PetscScalar));
304: } else {
305: PetscReal *cB = &cellBr[nc * offset];
307: DMFieldDSdot(cB, T->T[0], elem, nq, feDim, nc, PetscRealPart);
308: }
309: }
310: if (D) {
311: if (datatype == PETSC_SCALAR) {
312: PetscScalar *cD = &cellDs[nc * dim * offset];
314: DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), (PetscScalar));
315: for (p = 0; p < nq; p++) {
316: for (g = 0; g < nc; g++) {
317: PetscScalar vs[3];
319: for (d = 0; d < dimR; d++) {
320: vs[d] = 0.;
321: for (e = 0; e < dimR; e++) vs[d] += invJ[dimR * dimR * p + e * dimR + d] * cD[(nc * p + g) * dimR + e];
322: }
323: for (d = 0; d < dimR; d++) cD[(nc * p + g) * dimR + d] = vs[d];
324: }
325: }
326: } else {
327: PetscReal *cD = &cellDr[nc * dim * offset];
329: DMFieldDSdot(cD, T->T[1], elem, nq, feDim, (nc * dim), PetscRealPart);
330: for (p = 0; p < nq; p++) {
331: for (g = 0; g < nc; g++) {
332: for (d = 0; d < dimR; d++) {
333: v[d] = 0.;
334: for (e = 0; e < dimR; e++) v[d] += invJ[dimR * dimR * p + e * dimR + d] * cD[(nc * p + g) * dimR + e];
335: }
336: for (d = 0; d < dimR; d++) cD[(nc * p + g) * dimR + d] = v[d];
337: }
338: }
339: }
340: }
341: if (H) {
342: if (datatype == PETSC_SCALAR) {
343: PetscScalar *cH = &cellHs[nc * dim * dim * offset];
345: DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), (PetscScalar));
346: for (p = 0; p < nq; p++) {
347: for (g = 0; g < nc * dimR; g++) {
348: PetscScalar vs[3];
350: for (d = 0; d < dimR; d++) {
351: vs[d] = 0.;
352: for (e = 0; e < dimR; e++) vs[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[(nc * dimR * p + g) * dimR + e];
353: }
354: for (d = 0; d < dimR; d++) cH[(nc * dimR * p + g) * dimR + d] = vs[d];
355: }
356: for (g = 0; g < nc; g++) {
357: for (f = 0; f < dimR; f++) {
358: PetscScalar vs[3];
360: for (d = 0; d < dimR; d++) {
361: vs[d] = 0.;
362: for (e = 0; e < dimR; e++) vs[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[((nc * p + g) * dimR + e) * dimR + f];
363: }
364: for (d = 0; d < dimR; d++) cH[((nc * p + g) * dimR + d) * dimR + f] = vs[d];
365: }
366: }
367: }
368: } else {
369: PetscReal *cH = &cellHr[nc * dim * dim * offset];
371: DMFieldDSdot(cH, T->T[2], elem, nq, feDim, (nc * dim * dim), PetscRealPart);
372: for (p = 0; p < nq; p++) {
373: for (g = 0; g < nc * dimR; g++) {
374: for (d = 0; d < dimR; d++) {
375: v[d] = 0.;
376: for (e = 0; e < dimR; e++) v[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[(nc * dimR * p + g) * dimR + e];
377: }
378: for (d = 0; d < dimR; d++) cH[(nc * dimR * p + g) * dimR + d] = v[d];
379: }
380: for (g = 0; g < nc; g++) {
381: for (f = 0; f < dimR; f++) {
382: for (d = 0; d < dimR; d++) {
383: v[d] = 0.;
384: for (e = 0; e < dimR; e++) v[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[((nc * p + g) * dimR + e) * dimR + f];
385: }
386: for (d = 0; d < dimR; d++) cH[((nc * p + g) * dimR + d) * dimR + f] = v[d];
387: }
388: }
389: }
390: }
391: }
392: DMFieldRestoreClosure_Internal(field, c, &isDG, &closureSize, &array, &elem);
393: PetscTabulationDestroy(&T);
394: }
395: offset += nq;
396: }
397: {
398: MPI_Datatype origtype;
400: if (datatype == PETSC_SCALAR) {
401: origtype = MPIU_SCALAR;
402: } else {
403: origtype = MPIU_REAL;
404: }
405: if (B) {
406: MPI_Datatype Btype;
408: MPI_Type_contiguous(nc, origtype, &Btype);
409: MPI_Type_commit(&Btype);
410: PetscSFScatterBegin(cellSF, Btype, (datatype == PETSC_SCALAR) ? (void *)cellBs : (void *)cellBr, B);
411: PetscSFScatterEnd(cellSF, Btype, (datatype == PETSC_SCALAR) ? (void *)cellBs : (void *)cellBr, B);
412: MPI_Type_free(&Btype);
413: }
414: if (D) {
415: MPI_Datatype Dtype;
417: MPI_Type_contiguous(nc * dim, origtype, &Dtype);
418: MPI_Type_commit(&Dtype);
419: PetscSFScatterBegin(cellSF, Dtype, (datatype == PETSC_SCALAR) ? (void *)cellDs : (void *)cellDr, D);
420: PetscSFScatterEnd(cellSF, Dtype, (datatype == PETSC_SCALAR) ? (void *)cellDs : (void *)cellDr, D);
421: MPI_Type_free(&Dtype);
422: }
423: if (H) {
424: MPI_Datatype Htype;
426: MPI_Type_contiguous(nc * dim * dim, origtype, &Htype);
427: MPI_Type_commit(&Htype);
428: PetscSFScatterBegin(cellSF, Htype, (datatype == PETSC_SCALAR) ? (void *)cellHs : (void *)cellHr, H);
429: PetscSFScatterEnd(cellSF, Htype, (datatype == PETSC_SCALAR) ? (void *)cellHs : (void *)cellHr, H);
430: MPI_Type_free(&Htype);
431: }
432: }
433: PetscFree4(v, J, invJ, detJ);
434: PetscFree3(cellBr, cellDr, cellHr);
435: PetscFree3(cellBs, cellDs, cellHs);
436: PetscFree3(cellPoints, coordsReal, coordsRef);
437: MPI_Type_free(&pointType);
438: PetscSFDestroy(&cellSF);
439: return 0;
440: }
442: static PetscErrorCode DMFieldEvaluateFV_DS(DMField field, IS pointIS, PetscDataType type, void *B, void *D, void *H)
443: {
444: DMField_DS *dsfield = (DMField_DS *)field->data;
445: PetscInt h, imin;
446: PetscInt dim;
447: PetscClassId id;
448: PetscQuadrature quad = NULL;
449: PetscInt maxDegree;
450: PetscFEGeom *geom;
451: PetscInt Nq, Nc, dimC, qNc, N;
452: PetscInt numPoints;
453: void *qB = NULL, *qD = NULL, *qH = NULL;
454: const PetscReal *weights;
455: MPI_Datatype mpitype = type == PETSC_SCALAR ? MPIU_SCALAR : MPIU_REAL;
456: PetscObject disc;
457: DMField coordField;
459: Nc = field->numComponents;
460: DMGetCoordinateDim(field->dm, &dimC);
461: DMGetDimension(field->dm, &dim);
462: ISGetLocalSize(pointIS, &numPoints);
463: ISGetMinMax(pointIS, &imin, NULL);
464: for (h = 0; h < dsfield->height; h++) {
465: PetscInt hEnd;
467: DMPlexGetHeightStratum(field->dm, h, NULL, &hEnd);
468: if (imin < hEnd) break;
469: }
470: dim -= h;
471: DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc);
472: PetscObjectGetClassId(disc, &id);
474: DMGetCoordinateField(field->dm, &coordField);
475: DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree);
476: if (maxDegree <= 1) DMFieldCreateDefaultQuadrature(coordField, pointIS, &quad);
477: if (!quad) DMFieldCreateDefaultQuadrature(field, pointIS, &quad);
479: DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FALSE, &geom);
480: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &weights);
482: N = numPoints * Nq * Nc;
483: if (B) DMGetWorkArray(field->dm, N, mpitype, &qB);
484: if (D) DMGetWorkArray(field->dm, N * dimC, mpitype, &qD);
485: if (H) DMGetWorkArray(field->dm, N * dimC * dimC, mpitype, &qH);
486: DMFieldEvaluateFE(field, pointIS, quad, type, qB, qD, qH);
487: if (B) {
488: PetscInt i, j, k;
490: if (type == PETSC_SCALAR) {
491: PetscScalar *sB = (PetscScalar *)B;
492: PetscScalar *sqB = (PetscScalar *)qB;
494: for (i = 0; i < numPoints; i++) {
495: PetscReal vol = 0.;
497: for (j = 0; j < Nc; j++) sB[i * Nc + j] = 0.;
498: for (k = 0; k < Nq; k++) {
499: vol += geom->detJ[i * Nq + k] * weights[k];
500: for (j = 0; j < Nc; j++) sB[i * Nc + j] += geom->detJ[i * Nq + k] * weights[k] * sqB[(i * Nq + k) * Nc + j];
501: }
502: for (k = 0; k < Nc; k++) sB[i * Nc + k] /= vol;
503: }
504: } else {
505: PetscReal *rB = (PetscReal *)B;
506: PetscReal *rqB = (PetscReal *)qB;
508: for (i = 0; i < numPoints; i++) {
509: PetscReal vol = 0.;
511: for (j = 0; j < Nc; j++) rB[i * Nc + j] = 0.;
512: for (k = 0; k < Nq; k++) {
513: vol += geom->detJ[i * Nq + k] * weights[k];
514: for (j = 0; j < Nc; j++) rB[i * Nc + j] += weights[k] * rqB[(i * Nq + k) * Nc + j];
515: }
516: for (k = 0; k < Nc; k++) rB[i * Nc + k] /= vol;
517: }
518: }
519: }
520: if (D) {
521: PetscInt i, j, k, l, m;
523: if (type == PETSC_SCALAR) {
524: PetscScalar *sD = (PetscScalar *)D;
525: PetscScalar *sqD = (PetscScalar *)qD;
527: for (i = 0; i < numPoints; i++) {
528: PetscReal vol = 0.;
530: for (j = 0; j < Nc * dimC; j++) sD[i * Nc * dimC + j] = 0.;
531: for (k = 0; k < Nq; k++) {
532: vol += geom->detJ[i * Nq + k] * weights[k];
533: for (j = 0; j < Nc; j++) {
534: PetscScalar pD[3] = {0., 0., 0.};
536: for (l = 0; l < dimC; l++) {
537: for (m = 0; m < dim; m++) pD[l] += geom->invJ[((i * Nq + k) * dimC + m) * dimC + l] * sqD[((i * Nq + k) * Nc + j) * dim + m];
538: }
539: for (l = 0; l < dimC; l++) sD[(i * Nc + j) * dimC + l] += geom->detJ[i * Nq + k] * weights[k] * pD[l];
540: }
541: }
542: for (k = 0; k < Nc * dimC; k++) sD[i * Nc * dimC + k] /= vol;
543: }
544: } else {
545: PetscReal *rD = (PetscReal *)D;
546: PetscReal *rqD = (PetscReal *)qD;
548: for (i = 0; i < numPoints; i++) {
549: PetscReal vol = 0.;
551: for (j = 0; j < Nc * dimC; j++) rD[i * Nc * dimC + j] = 0.;
552: for (k = 0; k < Nq; k++) {
553: vol += geom->detJ[i * Nq + k] * weights[k];
554: for (j = 0; j < Nc; j++) {
555: PetscReal pD[3] = {0., 0., 0.};
557: for (l = 0; l < dimC; l++) {
558: for (m = 0; m < dim; m++) pD[l] += geom->invJ[((i * Nq + k) * dimC + m) * dimC + l] * rqD[((i * Nq + k) * Nc + j) * dim + m];
559: }
560: for (l = 0; l < dimC; l++) rD[(i * Nc + j) * dimC + l] += geom->detJ[i * Nq + k] * weights[k] * pD[l];
561: }
562: }
563: for (k = 0; k < Nc * dimC; k++) rD[i * Nc * dimC + k] /= vol;
564: }
565: }
566: }
567: if (H) {
568: PetscInt i, j, k, l, m, q, r;
570: if (type == PETSC_SCALAR) {
571: PetscScalar *sH = (PetscScalar *)H;
572: PetscScalar *sqH = (PetscScalar *)qH;
574: for (i = 0; i < numPoints; i++) {
575: PetscReal vol = 0.;
577: for (j = 0; j < Nc * dimC * dimC; j++) sH[i * Nc * dimC * dimC + j] = 0.;
578: for (k = 0; k < Nq; k++) {
579: const PetscReal *invJ = &geom->invJ[(i * Nq + k) * dimC * dimC];
581: vol += geom->detJ[i * Nq + k] * weights[k];
582: for (j = 0; j < Nc; j++) {
583: PetscScalar pH[3][3] = {
584: {0., 0., 0.},
585: {0., 0., 0.},
586: {0., 0., 0.}
587: };
588: const PetscScalar *spH = &sqH[((i * Nq + k) * Nc + j) * dimC * dimC];
590: for (l = 0; l < dimC; l++) {
591: for (m = 0; m < dimC; m++) {
592: for (q = 0; q < dim; q++) {
593: for (r = 0; r < dim; r++) pH[l][m] += invJ[q * dimC + l] * invJ[r * dimC + m] * spH[q * dim + r];
594: }
595: }
596: }
597: for (l = 0; l < dimC; l++) {
598: for (m = 0; m < dimC; m++) sH[(i * Nc + j) * dimC * dimC + l * dimC + m] += geom->detJ[i * Nq + k] * weights[k] * pH[l][m];
599: }
600: }
601: }
602: for (k = 0; k < Nc * dimC * dimC; k++) sH[i * Nc * dimC * dimC + k] /= vol;
603: }
604: } else {
605: PetscReal *rH = (PetscReal *)H;
606: PetscReal *rqH = (PetscReal *)qH;
608: for (i = 0; i < numPoints; i++) {
609: PetscReal vol = 0.;
611: for (j = 0; j < Nc * dimC * dimC; j++) rH[i * Nc * dimC * dimC + j] = 0.;
612: for (k = 0; k < Nq; k++) {
613: const PetscReal *invJ = &geom->invJ[(i * Nq + k) * dimC * dimC];
615: vol += geom->detJ[i * Nq + k] * weights[k];
616: for (j = 0; j < Nc; j++) {
617: PetscReal pH[3][3] = {
618: {0., 0., 0.},
619: {0., 0., 0.},
620: {0., 0., 0.}
621: };
622: const PetscReal *rpH = &rqH[((i * Nq + k) * Nc + j) * dimC * dimC];
624: for (l = 0; l < dimC; l++) {
625: for (m = 0; m < dimC; m++) {
626: for (q = 0; q < dim; q++) {
627: for (r = 0; r < dim; r++) pH[l][m] += invJ[q * dimC + l] * invJ[r * dimC + m] * rpH[q * dim + r];
628: }
629: }
630: }
631: for (l = 0; l < dimC; l++) {
632: for (m = 0; m < dimC; m++) rH[(i * Nc + j) * dimC * dimC + l * dimC + m] += geom->detJ[i * Nq + k] * weights[k] * pH[l][m];
633: }
634: }
635: }
636: for (k = 0; k < Nc * dimC * dimC; k++) rH[i * Nc * dimC * dimC + k] /= vol;
637: }
638: }
639: }
640: if (B) DMRestoreWorkArray(field->dm, N, mpitype, &qB);
641: if (D) DMRestoreWorkArray(field->dm, N * dimC, mpitype, &qD);
642: if (H) DMRestoreWorkArray(field->dm, N * dimC * dimC, mpitype, &qH);
643: PetscFEGeomDestroy(&geom);
644: PetscQuadratureDestroy(&quad);
645: return 0;
646: }
648: static PetscErrorCode DMFieldGetDegree_DS(DMField field, IS pointIS, PetscInt *minDegree, PetscInt *maxDegree)
649: {
650: DMField_DS *dsfield;
651: PetscObject disc;
652: PetscInt h, imin, imax;
653: PetscClassId id;
655: dsfield = (DMField_DS *)field->data;
656: ISGetMinMax(pointIS, &imin, &imax);
657: if (imin >= imax) {
658: h = 0;
659: } else {
660: for (h = 0; h < dsfield->height; h++) {
661: PetscInt hEnd;
663: DMPlexGetHeightStratum(field->dm, h, NULL, &hEnd);
664: if (imin < hEnd) break;
665: }
666: }
667: DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc);
668: PetscObjectGetClassId(disc, &id);
669: if (id == PETSCFE_CLASSID) {
670: PetscFE fe = (PetscFE)disc;
671: PetscSpace sp;
673: PetscFEGetBasisSpace(fe, &sp);
674: PetscSpaceGetDegree(sp, minDegree, maxDegree);
675: }
676: return 0;
677: }
679: PetscErrorCode DMFieldGetFVQuadrature_Internal(DMField field, IS pointIS, PetscQuadrature *quad)
680: {
681: DM dm = field->dm;
682: const PetscInt *points;
683: DMPolytopeType ct;
684: PetscInt dim, n;
685: PetscBool isplex;
687: PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isplex);
688: ISGetLocalSize(pointIS, &n);
689: if (isplex && n) {
690: DMGetDimension(dm, &dim);
691: ISGetIndices(pointIS, &points);
692: DMPlexGetCellType(dm, points[0], &ct);
693: switch (ct) {
694: case DM_POLYTOPE_TRIANGLE:
695: case DM_POLYTOPE_TETRAHEDRON:
696: PetscDTStroudConicalQuadrature(dim, 1, 1, -1.0, 1.0, quad);
697: break;
698: default:
699: PetscDTGaussTensorQuadrature(dim, 1, 1, -1.0, 1.0, quad);
700: }
701: ISRestoreIndices(pointIS, &points);
702: } else DMFieldCreateDefaultQuadrature(field, pointIS, quad);
703: return 0;
704: }
706: static PetscErrorCode DMFieldCreateDefaultQuadrature_DS(DMField field, IS pointIS, PetscQuadrature *quad)
707: {
708: PetscInt h, dim, imax, imin, cellHeight;
709: DM dm;
710: DMField_DS *dsfield;
711: PetscObject disc;
712: PetscFE fe;
713: PetscClassId id;
715: dm = field->dm;
716: dsfield = (DMField_DS *)field->data;
717: ISGetMinMax(pointIS, &imin, &imax);
718: DMGetDimension(dm, &dim);
719: for (h = 0; h <= dim; h++) {
720: PetscInt hStart, hEnd;
722: DMPlexGetHeightStratum(dm, h, &hStart, &hEnd);
723: if (imax >= hStart && imin < hEnd) break;
724: }
725: DMPlexGetVTKCellHeight(dm, &cellHeight);
726: h -= cellHeight;
727: *quad = NULL;
728: if (h < dsfield->height) {
729: DMFieldDSGetHeightDisc(field, h, dsfield->disc, &disc);
730: PetscObjectGetClassId(disc, &id);
731: if (id != PETSCFE_CLASSID) return 0;
732: fe = (PetscFE)disc;
733: PetscFEGetQuadrature(fe, quad);
734: PetscObjectReference((PetscObject)*quad);
735: }
736: return 0;
737: }
739: static PetscErrorCode DMFieldComputeFaceData_DS(DMField field, IS pointIS, PetscQuadrature quad, PetscFEGeom *geom)
740: {
741: const PetscInt *points;
742: PetscInt p, dim, dE, numFaces, Nq;
743: PetscInt maxDegree;
744: DMLabel depthLabel;
745: IS cellIS;
746: DM dm = field->dm;
748: dim = geom->dim;
749: dE = geom->dimEmbed;
750: DMPlexGetDepthLabel(dm, &depthLabel);
751: DMLabelGetStratumIS(depthLabel, dim + 1, &cellIS);
752: DMFieldGetDegree(field, cellIS, NULL, &maxDegree);
753: ISGetIndices(pointIS, &points);
754: numFaces = geom->numCells;
755: Nq = geom->numPoints;
756: /* First, set local faces and flip normals so that they are outward for the first supporting cell */
757: for (p = 0; p < numFaces; p++) {
758: PetscInt point = points[p];
759: PetscInt suppSize, s, coneSize, c, numChildren;
760: const PetscInt *supp, *cone, *ornt;
762: DMPlexGetTreeChildren(dm, point, &numChildren, NULL);
764: DMPlexGetSupportSize(dm, point, &suppSize);
766: if (!suppSize) continue;
767: DMPlexGetSupport(dm, point, &supp);
768: for (s = 0; s < suppSize; ++s) {
769: DMPlexGetConeSize(dm, supp[s], &coneSize);
770: DMPlexGetCone(dm, supp[s], &cone);
771: for (c = 0; c < coneSize; ++c)
772: if (cone[c] == point) break;
774: geom->face[p][s] = c;
775: }
776: DMPlexGetConeOrientation(dm, supp[0], &ornt);
777: if (ornt[geom->face[p][0]] < 0) {
778: PetscInt Np = geom->numPoints, q, dE = geom->dimEmbed, d;
780: for (q = 0; q < Np; ++q)
781: for (d = 0; d < dE; ++d) geom->n[(p * Np + q) * dE + d] = -geom->n[(p * Np + q) * dE + d];
782: }
783: }
784: if (maxDegree <= 1) {
785: PetscInt numCells, offset, *cells;
786: PetscFEGeom *cellGeom;
787: IS suppIS;
789: for (p = 0, numCells = 0; p < numFaces; p++) {
790: PetscInt point = points[p];
791: PetscInt numSupp, numChildren;
793: DMPlexGetTreeChildren(dm, point, &numChildren, NULL);
795: DMPlexGetSupportSize(dm, point, &numSupp);
797: numCells += numSupp;
798: }
799: PetscMalloc1(numCells, &cells);
800: for (p = 0, offset = 0; p < numFaces; p++) {
801: PetscInt point = points[p];
802: PetscInt numSupp, s;
803: const PetscInt *supp;
805: DMPlexGetSupportSize(dm, point, &numSupp);
806: DMPlexGetSupport(dm, point, &supp);
807: for (s = 0; s < numSupp; s++, offset++) cells[offset] = supp[s];
808: }
809: ISCreateGeneral(PETSC_COMM_SELF, numCells, cells, PETSC_USE_POINTER, &suppIS);
810: DMFieldCreateFEGeom(field, suppIS, quad, PETSC_FALSE, &cellGeom);
811: for (p = 0, offset = 0; p < numFaces; p++) {
812: PetscInt point = points[p];
813: PetscInt numSupp, s, q;
814: const PetscInt *supp;
816: DMPlexGetSupportSize(dm, point, &numSupp);
817: DMPlexGetSupport(dm, point, &supp);
818: for (s = 0; s < numSupp; s++, offset++) {
819: for (q = 0; q < Nq * dE * dE; q++) {
820: geom->suppJ[s][p * Nq * dE * dE + q] = cellGeom->J[offset * Nq * dE * dE + q];
821: geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q];
822: }
823: for (q = 0; q < Nq; q++) geom->suppDetJ[s][p * Nq + q] = cellGeom->detJ[offset * Nq + q];
824: }
825: }
826: PetscFEGeomDestroy(&cellGeom);
827: ISDestroy(&suppIS);
828: PetscFree(cells);
829: } else {
830: DMField_DS *dsfield = (DMField_DS *)field->data;
831: PetscObject faceDisc, cellDisc;
832: PetscClassId faceId, cellId;
833: PetscDualSpace dsp;
834: DM K;
835: DMPolytopeType ct;
836: PetscInt(*co)[2][3];
837: PetscInt coneSize;
838: PetscInt **counts;
839: PetscInt f, i, o, q, s;
840: const PetscInt *coneK;
841: PetscInt eStart, minOrient, maxOrient, numOrient;
842: PetscInt *orients;
843: PetscReal **orientPoints;
844: PetscReal *cellPoints;
845: PetscReal *dummyWeights;
846: PetscQuadrature cellQuad = NULL;
848: DMFieldDSGetHeightDisc(field, 1, dsfield->disc, &faceDisc);
849: DMFieldDSGetHeightDisc(field, 0, dsfield->disc, &cellDisc);
850: PetscObjectGetClassId(faceDisc, &faceId);
851: PetscObjectGetClassId(cellDisc, &cellId);
853: PetscFEGetDualSpace((PetscFE)cellDisc, &dsp);
854: PetscDualSpaceGetDM(dsp, &K);
855: DMPlexGetHeightStratum(K, 1, &eStart, NULL);
856: DMPlexGetCellType(K, eStart, &ct);
857: DMPlexGetConeSize(K, 0, &coneSize);
858: DMPlexGetCone(K, 0, &coneK);
859: PetscMalloc2(numFaces, &co, coneSize, &counts);
860: PetscMalloc1(dE * Nq, &cellPoints);
861: PetscMalloc1(Nq, &dummyWeights);
862: PetscQuadratureCreate(PETSC_COMM_SELF, &cellQuad);
863: PetscQuadratureSetData(cellQuad, dE, 1, Nq, cellPoints, dummyWeights);
864: minOrient = PETSC_MAX_INT;
865: maxOrient = PETSC_MIN_INT;
866: for (p = 0; p < numFaces; p++) { /* record the orientation of the facet wrt the support cells */
867: PetscInt point = points[p];
868: PetscInt numSupp, numChildren;
869: const PetscInt *supp;
871: DMPlexGetTreeChildren(dm, point, &numChildren, NULL);
873: DMPlexGetSupportSize(dm, point, &numSupp);
875: DMPlexGetSupport(dm, point, &supp);
876: for (s = 0; s < numSupp; s++) {
877: PetscInt cell = supp[s];
878: PetscInt numCone;
879: const PetscInt *cone, *orient;
881: DMPlexGetConeSize(dm, cell, &numCone);
883: DMPlexGetCone(dm, cell, &cone);
884: DMPlexGetConeOrientation(dm, cell, &orient);
885: for (f = 0; f < coneSize; f++) {
886: if (cone[f] == point) break;
887: }
888: co[p][s][0] = f;
889: co[p][s][1] = orient[f];
890: co[p][s][2] = cell;
891: minOrient = PetscMin(minOrient, orient[f]);
892: maxOrient = PetscMax(maxOrient, orient[f]);
893: }
894: for (; s < 2; s++) {
895: co[p][s][0] = -1;
896: co[p][s][1] = -1;
897: co[p][s][2] = -1;
898: }
899: }
900: numOrient = maxOrient + 1 - minOrient;
901: DMPlexGetCone(K, 0, &coneK);
902: /* count all (face,orientation) doubles that appear */
903: PetscCalloc2(numOrient, &orients, numOrient, &orientPoints);
904: for (f = 0; f < coneSize; f++) PetscCalloc1(numOrient + 1, &counts[f]);
905: for (p = 0; p < numFaces; p++) {
906: for (s = 0; s < 2; s++) {
907: if (co[p][s][0] >= 0) {
908: counts[co[p][s][0]][co[p][s][1] - minOrient]++;
909: orients[co[p][s][1] - minOrient]++;
910: }
911: }
912: }
913: for (o = 0; o < numOrient; o++) {
914: if (orients[o]) {
915: PetscInt orient = o + minOrient;
916: PetscInt q;
918: PetscMalloc1(Nq * dim, &orientPoints[o]);
919: /* rotate the quadrature points appropriately */
920: switch (ct) {
921: case DM_POLYTOPE_POINT:
922: break;
923: case DM_POLYTOPE_SEGMENT:
924: if (orient == -2 || orient == 1) {
925: for (q = 0; q < Nq; q++) orientPoints[o][q] = -geom->xi[q];
926: } else {
927: for (q = 0; q < Nq; q++) orientPoints[o][q] = geom->xi[q];
928: }
929: break;
930: case DM_POLYTOPE_TRIANGLE:
931: for (q = 0; q < Nq; q++) {
932: PetscReal lambda[3];
933: PetscReal lambdao[3];
935: /* convert to barycentric */
936: lambda[0] = -(geom->xi[2 * q] + geom->xi[2 * q + 1]) / 2.;
937: lambda[1] = (geom->xi[2 * q] + 1.) / 2.;
938: lambda[2] = (geom->xi[2 * q + 1] + 1.) / 2.;
939: if (orient >= 0) {
940: for (i = 0; i < 3; i++) lambdao[i] = lambda[(orient + i) % 3];
941: } else {
942: for (i = 0; i < 3; i++) lambdao[i] = lambda[(-(orient + i) + 3) % 3];
943: }
944: /* convert to coordinates */
945: orientPoints[o][2 * q + 0] = -(lambdao[0] + lambdao[2]) + lambdao[1];
946: orientPoints[o][2 * q + 1] = -(lambdao[0] + lambdao[1]) + lambdao[2];
947: }
948: break;
949: case DM_POLYTOPE_QUADRILATERAL:
950: for (q = 0; q < Nq; q++) {
951: PetscReal xi[2], xio[2];
952: PetscInt oabs = (orient >= 0) ? orient : -(orient + 1);
954: xi[0] = geom->xi[2 * q];
955: xi[1] = geom->xi[2 * q + 1];
956: switch (oabs) {
957: case 1:
958: xio[0] = xi[1];
959: xio[1] = -xi[0];
960: break;
961: case 2:
962: xio[0] = -xi[0];
963: xio[1] = -xi[1];
964: break;
965: case 3:
966: xio[0] = -xi[1];
967: xio[1] = xi[0];
968: break;
969: case 0:
970: default:
971: xio[0] = xi[0];
972: xio[1] = xi[1];
973: break;
974: }
975: if (orient < 0) xio[0] = -xio[0];
976: orientPoints[o][2 * q + 0] = xio[0];
977: orientPoints[o][2 * q + 1] = xio[1];
978: }
979: break;
980: default:
981: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s not yet supported", DMPolytopeTypes[ct]);
982: }
983: }
984: }
985: for (f = 0; f < coneSize; f++) {
986: PetscInt face = coneK[f];
987: PetscReal v0[3];
988: PetscReal J[9], detJ;
989: PetscInt numCells, offset;
990: PetscInt *cells;
991: IS suppIS;
993: DMPlexComputeCellGeometryFEM(K, face, NULL, v0, J, NULL, &detJ);
994: for (o = 0; o <= numOrient; o++) {
995: PetscFEGeom *cellGeom;
997: if (!counts[f][o]) continue;
998: /* If this (face,orientation) double appears,
999: * convert the face quadrature points into volume quadrature points */
1000: for (q = 0; q < Nq; q++) {
1001: PetscReal xi0[3] = {-1., -1., -1.};
1003: CoordinatesRefToReal(dE, dim, xi0, v0, J, &orientPoints[o][dim * q + 0], &cellPoints[dE * q + 0]);
1004: }
1005: for (p = 0, numCells = 0; p < numFaces; p++) {
1006: for (s = 0; s < 2; s++) {
1007: if (co[p][s][0] == f && co[p][s][1] == o + minOrient) numCells++;
1008: }
1009: }
1010: PetscMalloc1(numCells, &cells);
1011: for (p = 0, offset = 0; p < numFaces; p++) {
1012: for (s = 0; s < 2; s++) {
1013: if (co[p][s][0] == f && co[p][s][1] == o + minOrient) cells[offset++] = co[p][s][2];
1014: }
1015: }
1016: ISCreateGeneral(PETSC_COMM_SELF, numCells, cells, PETSC_USE_POINTER, &suppIS);
1017: DMFieldCreateFEGeom(field, suppIS, cellQuad, PETSC_FALSE, &cellGeom);
1018: for (p = 0, offset = 0; p < numFaces; p++) {
1019: for (s = 0; s < 2; s++) {
1020: if (co[p][s][0] == f && co[p][s][1] == o + minOrient) {
1021: for (q = 0; q < Nq * dE * dE; q++) {
1022: geom->suppJ[s][p * Nq * dE * dE + q] = cellGeom->J[offset * Nq * dE * dE + q];
1023: geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q];
1024: }
1025: for (q = 0; q < Nq; q++) geom->suppDetJ[s][p * Nq + q] = cellGeom->detJ[offset * Nq + q];
1026: offset++;
1027: }
1028: }
1029: }
1030: PetscFEGeomDestroy(&cellGeom);
1031: ISDestroy(&suppIS);
1032: PetscFree(cells);
1033: }
1034: }
1035: for (o = 0; o < numOrient; o++) {
1036: if (orients[o]) PetscFree(orientPoints[o]);
1037: }
1038: PetscFree2(orients, orientPoints);
1039: PetscQuadratureDestroy(&cellQuad);
1040: for (f = 0; f < coneSize; f++) PetscFree(counts[f]);
1041: PetscFree2(co, counts);
1042: }
1043: ISRestoreIndices(pointIS, &points);
1044: ISDestroy(&cellIS);
1045: return 0;
1046: }
1048: static PetscErrorCode DMFieldInitialize_DS(DMField field)
1049: {
1050: field->ops->destroy = DMFieldDestroy_DS;
1051: field->ops->evaluate = DMFieldEvaluate_DS;
1052: field->ops->evaluateFE = DMFieldEvaluateFE_DS;
1053: field->ops->evaluateFV = DMFieldEvaluateFV_DS;
1054: field->ops->getDegree = DMFieldGetDegree_DS;
1055: field->ops->createDefaultQuadrature = DMFieldCreateDefaultQuadrature_DS;
1056: field->ops->view = DMFieldView_DS;
1057: field->ops->computeFaceData = DMFieldComputeFaceData_DS;
1058: return 0;
1059: }
1061: PETSC_INTERN PetscErrorCode DMFieldCreate_DS(DMField field)
1062: {
1063: DMField_DS *dsfield;
1065: PetscNew(&dsfield);
1066: field->data = dsfield;
1067: DMFieldInitialize_DS(field);
1068: return 0;
1069: }
1071: PetscErrorCode DMFieldCreateDSWithDG(DM dm, DM dmDG, PetscInt fieldNum, Vec vec, Vec vecDG, DMField *field)
1072: {
1073: DMField b;
1074: DMField_DS *dsfield;
1075: PetscObject disc = NULL, discDG = NULL;
1076: PetscSection section;
1077: PetscBool isContainer = PETSC_FALSE;
1078: PetscClassId id = -1;
1079: PetscInt numComponents = -1, dsNumFields;
1081: DMGetLocalSection(dm, §ion);
1082: PetscSectionGetFieldComponents(section, fieldNum, &numComponents);
1083: DMGetNumFields(dm, &dsNumFields);
1084: if (dsNumFields) DMGetField(dm, fieldNum, NULL, &disc);
1085: if (dsNumFields && dmDG) {
1086: DMGetField(dmDG, fieldNum, NULL, &discDG);
1087: PetscObjectReference(discDG);
1088: }
1089: if (disc) {
1090: PetscObjectGetClassId(disc, &id);
1091: isContainer = (id == PETSC_CONTAINER_CLASSID) ? PETSC_TRUE : PETSC_FALSE;
1092: }
1093: if (!disc || isContainer) {
1094: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
1095: PetscFE fe;
1096: DMPolytopeType ct, locct = DM_POLYTOPE_UNKNOWN;
1097: PetscInt dim, cStart, cEnd, cellHeight;
1099: DMPlexGetVTKCellHeight(dm, &cellHeight);
1100: DMGetDimension(dm, &dim);
1101: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
1102: if (cEnd > cStart) DMPlexGetCellType(dm, cStart, &locct);
1103: MPI_Allreduce(&locct, &ct, 1, MPI_INT, MPI_MIN, comm);
1104: PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, numComponents, ct, 1, PETSC_DETERMINE, &fe);
1105: PetscFEViewFromOptions(fe, NULL, "-field_fe_view");
1106: disc = (PetscObject)fe;
1107: } else PetscObjectReference(disc);
1108: PetscObjectGetClassId(disc, &id);
1109: if (id == PETSCFE_CLASSID) PetscFEGetNumComponents((PetscFE)disc, &numComponents);
1110: else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot determine number of discretization components");
1111: DMFieldCreate(dm, numComponents, DMFIELD_VERTEX, &b);
1112: DMFieldSetType(b, DMFIELDDS);
1113: dsfield = (DMField_DS *)b->data;
1114: dsfield->fieldNum = fieldNum;
1115: DMGetDimension(dm, &dsfield->height);
1116: dsfield->height++;
1117: PetscCalloc1(dsfield->height, &dsfield->disc);
1118: dsfield->disc[0] = disc;
1119: PetscObjectReference((PetscObject)vec);
1120: dsfield->vec = vec;
1121: if (dmDG) {
1122: dsfield->dmDG = dmDG;
1123: PetscCalloc1(dsfield->height, &dsfield->discDG);
1124: dsfield->discDG[0] = discDG;
1125: PetscObjectReference((PetscObject)vecDG);
1126: dsfield->vecDG = vecDG;
1127: }
1128: *field = b;
1129: return 0;
1130: }
1132: PetscErrorCode DMFieldCreateDS(DM dm, PetscInt fieldNum, Vec vec, DMField *field)
1133: {
1134: DMFieldCreateDSWithDG(dm, NULL, fieldNum, vec, NULL, field);
1135: return 0;
1136: }