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: {
9: PetscInt fieldNum;
10: Vec vec;
11: PetscInt height;
12: PetscObject *disc;
13: PetscBool multifieldVec;
14: }
15: DMField_DS;
17: static PetscErrorCode DMFieldDestroy_DS(DMField field)
18: {
19: DMField_DS *dsfield;
20: PetscInt i;
24: dsfield = (DMField_DS *) field->data;
25: VecDestroy(&dsfield->vec);
26: for (i = 0; i < dsfield->height; i++) {
27: PetscObjectDereference(dsfield->disc[i]);
28: }
29: PetscFree(dsfield->disc);
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: PetscBool iascii;
38: PetscObject disc;
42: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
43: disc = dsfield->disc[0];
44: if (iascii) {
45: PetscViewerASCIIPrintf(viewer, "PetscDS field %D\n",dsfield->fieldNum);
46: PetscViewerASCIIPushTab(viewer);
47: PetscObjectView(disc,viewer);
48: PetscViewerASCIIPopTab(viewer);
49: }
50: PetscViewerASCIIPushTab(viewer);
51: if (dsfield->multifieldVec) {
52: SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"View of subfield not implemented yet");
53: } else {
54: VecView(dsfield->vec,viewer);
55: }
56: PetscViewerASCIIPopTab(viewer);
57: return(0);
58: }
60: static PetscErrorCode DMFieldDSGetHeightDisc(DMField field, PetscInt height, PetscObject *disc)
61: {
62: DMField_DS *dsfield = (DMField_DS *) field->data;
66: if (!dsfield->disc[height]) {
67: PetscClassId id;
69: PetscObjectGetClassId(dsfield->disc[0],&id);
70: if (id == PETSCFE_CLASSID) {
71: PetscFE fe = (PetscFE) dsfield->disc[0];
73: PetscFECreateHeightTrace(fe,height,(PetscFE *)&dsfield->disc[height]);
74: }
75: }
76: *disc = dsfield->disc[height];
77: return(0);
78: }
80: /* y[m,c] = A[m,n,c] . b[n] */
81: #define DMFieldDSdot(y,A,b,m,n,c,cast) \
82: do { \
83: PetscInt _i, _j, _k; \
84: for (_i = 0; _i < (m); _i++) { \
85: for (_k = 0; _k < (c); _k++) { \
86: (y)[_i * (c) + _k] = 0.; \
87: } \
88: for (_j = 0; _j < (n); _j++) { \
89: for (_k = 0; _k < (c); _k++) { \
90: (y)[_i * (c) + _k] += (A)[(_i * (n) + _j) * (c) + _k] * cast((b)[_j]); \
91: } \
92: } \
93: } \
94: } while (0)
96: /* TODO: Reorganize interface so that I can reuse a tabulation rather than mallocing each time */
97: static PetscErrorCode DMFieldEvaluateFE_DS(DMField field, IS pointIS, PetscQuadrature quad, PetscDataType type, void *B, void *D, void *H)
98: {
99: DMField_DS *dsfield = (DMField_DS *) field->data;
100: DM dm;
101: PetscObject disc;
102: PetscClassId classid;
103: PetscInt nq, nc, dim, meshDim, numCells;
104: PetscSection section;
105: const PetscReal *qpoints;
106: PetscBool isStride;
107: const PetscInt *points = NULL;
108: PetscInt sfirst = -1, stride = -1;
109: PetscErrorCode ierr;
112: dm = field->dm;
113: nc = field->numComponents;
114: PetscQuadratureGetData(quad,&dim,NULL,&nq,&qpoints,NULL);
115: DMFieldDSGetHeightDisc(field,dsfield->height - 1 - dim,&disc);
116: DMGetDimension(dm,&meshDim);
117: DMGetLocalSection(dm,§ion);
118: PetscSectionGetField(section,dsfield->fieldNum,§ion);
119: PetscObjectGetClassId(disc,&classid);
120: /* TODO: batch */
121: PetscObjectTypeCompare((PetscObject)pointIS,ISSTRIDE,&isStride);
122: ISGetLocalSize(pointIS,&numCells);
123: if (isStride) {
124: ISStrideGetInfo(pointIS,&sfirst,&stride);
125: } else {
126: ISGetIndices(pointIS,&points);
127: }
128: if (classid == PETSCFE_CLASSID) {
129: PetscFE fe = (PetscFE) disc;
130: PetscInt feDim, i;
131: PetscInt K = H ? 2 : (D ? 1 : (B ? 0 : -1));
132: PetscTabulation T;
134: PetscFEGetDimension(fe,&feDim);
135: PetscFECreateTabulation(fe,1,nq,qpoints,K,&T);
136: for (i = 0; i < numCells; i++) {
137: PetscInt c = isStride ? (sfirst + i * stride) : points[i];
138: PetscInt closureSize;
139: PetscScalar *elem = NULL;
141: DMPlexVecGetClosure(dm,section,dsfield->vec,c,&closureSize,&elem);
142: if (B) {
143: /* field[c] = T[q,b,c] . coef[b], so v[c] = T[q,b,c] . coords[b] */
144: if (type == PETSC_SCALAR) {
145: PetscScalar *cB = &((PetscScalar *) B)[nc * nq * i];
147: DMFieldDSdot(cB,T->T[0],elem,nq,feDim,nc,(PetscScalar));
148: } else {
149: PetscReal *cB = &((PetscReal *) B)[nc * nq * i];
151: DMFieldDSdot(cB,T->T[0],elem,nq,feDim,nc,PetscRealPart);
152: }
153: }
154: if (D) {
155: if (type == PETSC_SCALAR) {
156: PetscScalar *cD = &((PetscScalar *) D)[nc * nq * dim * i];
158: DMFieldDSdot(cD,T->T[1],elem,nq,feDim,(nc * dim),(PetscScalar));
159: } else {
160: PetscReal *cD = &((PetscReal *) D)[nc * nq * dim * i];
162: DMFieldDSdot(cD,T->T[1],elem,nq,feDim,(nc * dim),PetscRealPart);
163: }
164: }
165: if (H) {
166: if (type == PETSC_SCALAR) {
167: PetscScalar *cH = &((PetscScalar *) H)[nc * nq * dim * dim * i];
169: DMFieldDSdot(cH,T->T[2],elem,nq,feDim,(nc * dim * dim),(PetscScalar));
170: } else {
171: PetscReal *cH = &((PetscReal *) H)[nc * nq * dim * dim * i];
173: DMFieldDSdot(cH,T->T[2],elem,nq,feDim,(nc * dim * dim),PetscRealPart);
174: }
175: }
176: DMPlexVecRestoreClosure(dm,section,dsfield->vec,c,&closureSize,&elem);
177: }
178: PetscTabulationDestroy(&T);
179: } else SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented");
180: if (!isStride) {
181: ISRestoreIndices(pointIS,&points);
182: }
183: return(0);
184: }
186: static PetscErrorCode DMFieldEvaluate_DS(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
187: {
188: DMField_DS *dsfield = (DMField_DS *) field->data;
189: PetscSF cellSF = NULL;
190: const PetscSFNode *cells;
191: PetscInt c, nFound, numCells, feDim, nc;
192: const PetscInt *cellDegrees;
193: const PetscScalar *pointsArray;
194: PetscScalar *cellPoints;
195: PetscInt gatherSize, gatherMax;
196: PetscInt dim, dimR, offset;
197: MPI_Datatype pointType;
198: PetscObject cellDisc;
199: PetscFE cellFE;
200: PetscClassId discID;
201: PetscReal *coordsReal, *coordsRef;
202: PetscSection section;
203: PetscScalar *cellBs = NULL, *cellDs = NULL, *cellHs = NULL;
204: PetscReal *cellBr = NULL, *cellDr = NULL, *cellHr = NULL;
205: PetscReal *v, *J, *invJ, *detJ;
206: PetscErrorCode ierr;
209: nc = field->numComponents;
210: DMGetLocalSection(field->dm,§ion);
211: DMFieldDSGetHeightDisc(field,0,&cellDisc);
212: PetscObjectGetClassId(cellDisc, &discID);
213: if (discID != PETSCFE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Discretization type not supported\n");
214: cellFE = (PetscFE) cellDisc;
215: PetscFEGetDimension(cellFE,&feDim);
216: DMGetCoordinateDim(field->dm, &dim);
217: DMGetDimension(field->dm, &dimR);
218: DMLocatePoints(field->dm, points, DM_POINTLOCATION_NONE, &cellSF);
219: PetscSFGetGraph(cellSF, &numCells, &nFound, NULL, &cells);
220: for (c = 0; c < nFound; c++) {
221: if (cells[c].index < 0) SETERRQ1(PetscObjectComm((PetscObject)points),PETSC_ERR_ARG_WRONG, "Point %D could not be located\n", c);
222: }
223: PetscSFComputeDegreeBegin(cellSF,&cellDegrees);
224: PetscSFComputeDegreeEnd(cellSF,&cellDegrees);
225: for (c = 0, gatherSize = 0, gatherMax = 0; c < numCells; c++) {
226: gatherMax = PetscMax(gatherMax,cellDegrees[c]);
227: gatherSize += cellDegrees[c];
228: }
229: PetscMalloc3(gatherSize*dim,&cellPoints,gatherMax*dim,&coordsReal,gatherMax*dimR,&coordsRef);
230: PetscMalloc4(gatherMax*dimR,&v,gatherMax*dimR*dimR,&J,gatherMax*dimR*dimR,&invJ,gatherMax,&detJ);
231: if (datatype == PETSC_SCALAR) {
232: PetscMalloc3(B ? nc * gatherSize : 0, &cellBs, D ? nc * dim * gatherSize : 0, &cellDs, H ? nc * dim * dim * gatherSize : 0, &cellHs);
233: } else {
234: PetscMalloc3(B ? nc * gatherSize : 0, &cellBr, D ? nc * dim * gatherSize : 0, &cellDr, H ? nc * dim * dim * gatherSize : 0, &cellHr);
235: }
237: MPI_Type_contiguous(dim,MPIU_SCALAR,&pointType);
238: MPI_Type_commit(&pointType);
239: VecGetArrayRead(points,&pointsArray);
240: PetscSFGatherBegin(cellSF, pointType, pointsArray, cellPoints);
241: PetscSFGatherEnd(cellSF, pointType, pointsArray, cellPoints);
242: VecRestoreArrayRead(points,&pointsArray);
243: for (c = 0, offset = 0; c < numCells; c++) {
244: PetscInt nq = cellDegrees[c], p;
246: if (nq) {
247: PetscInt K = H ? 2 : (D ? 1 : (B ? 0 : -1));
248: PetscTabulation T;
249: PetscInt closureSize;
250: PetscScalar *elem = NULL;
251: PetscReal *quadPoints;
252: PetscQuadrature quad;
253: PetscInt d, e, f, g;
255: for (p = 0; p < dim * nq; p++) coordsReal[p] = PetscRealPart(cellPoints[dim * offset + p]);
256: DMPlexCoordinatesToReference(field->dm, c, nq, coordsReal, coordsRef);
257: PetscFECreateTabulation(cellFE,1,nq,coordsRef,K,&T);
258: PetscQuadratureCreate(PETSC_COMM_SELF, &quad);
259: PetscMalloc1(dimR * nq, &quadPoints);
260: for (p = 0; p < dimR * nq; p++) quadPoints[p] = coordsRef[p];
261: PetscQuadratureSetData(quad, dimR, 0, nq, quadPoints, NULL);
262: DMPlexComputeCellGeometryFEM(field->dm, c, quad, v, J, invJ, detJ);
263: PetscQuadratureDestroy(&quad);
264: DMPlexVecGetClosure(field->dm,section,dsfield->vec,c,&closureSize,&elem);
265: if (B) {
266: if (datatype == PETSC_SCALAR) {
267: PetscScalar *cB = &cellBs[nc * offset];
269: DMFieldDSdot(cB,T->T[0],elem,nq,feDim,nc,(PetscScalar));
270: } else {
271: PetscReal *cB = &cellBr[nc * offset];
273: DMFieldDSdot(cB,T->T[0],elem,nq,feDim,nc,PetscRealPart);
274: }
275: }
276: if (D) {
277: if (datatype == PETSC_SCALAR) {
278: PetscScalar *cD = &cellDs[nc * dim * offset];
280: DMFieldDSdot(cD,T->T[1],elem,nq,feDim,(nc * dim),(PetscScalar));
281: for (p = 0; p < nq; p++) {
282: for (g = 0; g < nc; g++) {
283: PetscScalar vs[3];
285: for (d = 0; d < dimR; d++) {
286: vs[d] = 0.;
287: for (e = 0; e < dimR; e++) {
288: vs[d] += invJ[dimR * dimR * p + e * dimR + d] * cD[(nc * p + g) * dimR + e];
289: }
290: }
291: for (d = 0; d < dimR; d++) {
292: cD[(nc * p + g) * dimR + d] = vs[d];
293: }
294: }
295: }
296: } else {
297: PetscReal *cD = &cellDr[nc * dim * offset];
299: DMFieldDSdot(cD,T->T[1],elem,nq,feDim,(nc * dim),PetscRealPart);
300: for (p = 0; p < nq; p++) {
301: for (g = 0; g < nc; g++) {
302: for (d = 0; d < dimR; d++) {
303: v[d] = 0.;
304: for (e = 0; e < dimR; e++) {
305: v[d] += invJ[dimR * dimR * p + e * dimR + d] * cD[(nc * p + g) * dimR + e];
306: }
307: }
308: for (d = 0; d < dimR; d++) {
309: cD[(nc * p + g) * dimR + d] = v[d];
310: }
311: }
312: }
313: }
314: }
315: if (H) {
316: if (datatype == PETSC_SCALAR) {
317: PetscScalar *cH = &cellHs[nc * dim * dim * offset];
319: DMFieldDSdot(cH,T->T[2],elem,nq,feDim,(nc * dim * dim),(PetscScalar));
320: for (p = 0; p < nq; p++) {
321: for (g = 0; g < nc * dimR; g++) {
322: PetscScalar vs[3];
324: for (d = 0; d < dimR; d++) {
325: vs[d] = 0.;
326: for (e = 0; e < dimR; e++) {
327: vs[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[(nc * dimR * p + g) * dimR + e];
328: }
329: }
330: for (d = 0; d < dimR; d++) {
331: cH[(nc * dimR * p + g) * dimR + d] = vs[d];
332: }
333: }
334: for (g = 0; g < nc; g++) {
335: for (f = 0; f < dimR; f++) {
336: PetscScalar vs[3];
338: for (d = 0; d < dimR; d++) {
339: vs[d] = 0.;
340: for (e = 0; e < dimR; e++) {
341: vs[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[((nc * p + g) * dimR + e) * dimR + f];
342: }
343: }
344: for (d = 0; d < dimR; d++) {
345: cH[((nc * p + g) * dimR + d) * dimR + f] = vs[d];
346: }
347: }
348: }
349: }
350: } else {
351: PetscReal *cH = &cellHr[nc * dim * dim * offset];
353: DMFieldDSdot(cH,T->T[2],elem,nq,feDim,(nc * dim * dim),PetscRealPart);
354: for (p = 0; p < nq; p++) {
355: for (g = 0; g < nc * dimR; g++) {
356: for (d = 0; d < dimR; d++) {
357: v[d] = 0.;
358: for (e = 0; e < dimR; e++) {
359: v[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[(nc * dimR * p + g) * dimR + e];
360: }
361: }
362: for (d = 0; d < dimR; d++) {
363: cH[(nc * dimR * p + g) * dimR + d] = v[d];
364: }
365: }
366: for (g = 0; g < nc; g++) {
367: for (f = 0; f < dimR; f++) {
368: for (d = 0; d < dimR; d++) {
369: v[d] = 0.;
370: for (e = 0; e < dimR; e++) {
371: v[d] += invJ[dimR * dimR * p + e * dimR + d] * cH[((nc * p + g) * dimR + e) * dimR + f];
372: }
373: }
374: for (d = 0; d < dimR; d++) {
375: cH[((nc * p + g) * dimR + d) * dimR + f] = v[d];
376: }
377: }
378: }
379: }
380: }
381: }
382: DMPlexVecRestoreClosure(field->dm,section,dsfield->vec,c,&closureSize,&elem);
383: PetscTabulationDestroy(&T);
384: }
385: offset += nq;
386: }
387: {
388: MPI_Datatype origtype;
390: if (datatype == PETSC_SCALAR) {
391: origtype = MPIU_SCALAR;
392: } else {
393: origtype = MPIU_REAL;
394: }
395: if (B) {
396: MPI_Datatype Btype;
398: MPI_Type_contiguous(nc, origtype, &Btype);
399: MPI_Type_commit(&Btype);
400: PetscSFScatterBegin(cellSF,Btype,(datatype == PETSC_SCALAR) ? (void *) cellBs : (void *) cellBr, B);
401: PetscSFScatterEnd(cellSF,Btype,(datatype == PETSC_SCALAR) ? (void *) cellBs : (void *) cellBr, B);
402: MPI_Type_free(&Btype);
403: }
404: if (D) {
405: MPI_Datatype Dtype;
407: MPI_Type_contiguous(nc * dim, origtype, &Dtype);
408: MPI_Type_commit(&Dtype);
409: PetscSFScatterBegin(cellSF,Dtype,(datatype == PETSC_SCALAR) ? (void *) cellDs : (void *) cellDr, D);
410: PetscSFScatterEnd(cellSF,Dtype,(datatype == PETSC_SCALAR) ? (void *) cellDs : (void *) cellDr, D);
411: MPI_Type_free(&Dtype);
412: }
413: if (H) {
414: MPI_Datatype Htype;
416: MPI_Type_contiguous(nc * dim * dim, origtype, &Htype);
417: MPI_Type_commit(&Htype);
418: PetscSFScatterBegin(cellSF,Htype,(datatype == PETSC_SCALAR) ? (void *) cellHs : (void *) cellHr, H);
419: PetscSFScatterEnd(cellSF,Htype,(datatype == PETSC_SCALAR) ? (void *) cellHs : (void *) cellHr, H);
420: MPI_Type_free(&Htype);
421: }
422: }
423: PetscFree4(v,J,invJ,detJ);
424: PetscFree3(cellBr, cellDr, cellHr);
425: PetscFree3(cellBs, cellDs, cellHs);
426: PetscFree3(cellPoints,coordsReal,coordsRef);
427: MPI_Type_free(&pointType);
428: PetscSFDestroy(&cellSF);
429: return(0);
430: }
432: static PetscErrorCode DMFieldEvaluateFV_DS(DMField field, IS pointIS, PetscDataType type, void *B, void *D, void *H)
433: {
434: DMField_DS *dsfield = (DMField_DS *) field->data;
435: PetscInt h, imin;
436: PetscInt dim;
437: PetscClassId id;
438: PetscQuadrature quad = NULL;
439: PetscInt maxDegree;
440: PetscFEGeom *geom;
441: PetscInt Nq, Nc, dimC, qNc, N;
442: PetscInt numPoints;
443: void *qB = NULL, *qD = NULL, *qH = NULL;
444: const PetscReal *weights;
445: MPI_Datatype mpitype = type == PETSC_SCALAR ? MPIU_SCALAR : MPIU_REAL;
446: PetscObject disc;
447: DMField coordField;
448: PetscErrorCode ierr;
451: Nc = field->numComponents;
452: DMGetCoordinateDim(field->dm, &dimC);
453: DMGetDimension(field->dm, &dim);
454: ISGetLocalSize(pointIS, &numPoints);
455: ISGetMinMax(pointIS,&imin,NULL);
456: for (h = 0; h < dsfield->height; h++) {
457: PetscInt hEnd;
459: DMPlexGetHeightStratum(field->dm,h,NULL,&hEnd);
460: if (imin < hEnd) break;
461: }
462: dim -= h;
463: DMFieldDSGetHeightDisc(field,h,&disc);
464: PetscObjectGetClassId(disc,&id);
465: if (id != PETSCFE_CLASSID) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Discretization not supported\n");
466: DMGetCoordinateField(field->dm, &coordField);
467: DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree);
468: if (maxDegree <= 1) {
469: DMFieldCreateDefaultQuadrature(coordField, pointIS, &quad);
470: }
471: if (!quad) {DMFieldCreateDefaultQuadrature(field, pointIS, &quad);}
472: if (!quad) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not determine quadrature for cell averages\n");
473: DMFieldCreateFEGeom(coordField,pointIS,quad,PETSC_FALSE,&geom);
474: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &weights);
475: if (qNc != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected scalar quadrature components\n");
476: N = numPoints * Nq * Nc;
477: if (B) {DMGetWorkArray(field->dm, N, mpitype, &qB);}
478: if (D) {DMGetWorkArray(field->dm, N * dimC, mpitype, &qD);}
479: if (H) {DMGetWorkArray(field->dm, N * dimC * dimC, mpitype, &qH);}
480: DMFieldEvaluateFE(field,pointIS,quad,type,qB,qD,qH);
481: if (B) {
482: PetscInt i, j, k;
484: if (type == PETSC_SCALAR) {
485: PetscScalar * sB = (PetscScalar *) B;
486: PetscScalar * sqB = (PetscScalar *) qB;
488: for (i = 0; i < numPoints; i++) {
489: PetscReal vol = 0.;
491: for (j = 0; j < Nc; j++) {sB[i * Nc + j] = 0.;}
492: for (k = 0; k < Nq; k++) {
493: vol += geom->detJ[i * Nq + k] * weights[k];
494: for (j = 0; j < Nc; j++) {
495: sB[i * Nc + j] += geom->detJ[i * Nq + k] * weights[k] * sqB[ (i * Nq + k) * Nc + j];
496: }
497: }
498: for (k = 0; k < Nc; k++) sB[i * Nc + k] /= vol;
499: }
500: } else {
501: PetscReal * rB = (PetscReal *) B;
502: PetscReal * rqB = (PetscReal *) qB;
504: for (i = 0; i < numPoints; i++) {
505: PetscReal vol = 0.;
507: for (j = 0; j < Nc; j++) {rB[i * Nc + j] = 0.;}
508: for (k = 0; k < Nq; k++) {
509: vol += geom->detJ[i * Nq + k] * weights[k];
510: for (j = 0; j < Nc; j++) {
511: rB[i * Nc + j] += weights[k] * rqB[ (i * Nq + k) * Nc + j];
512: }
513: }
514: for (k = 0; k < Nc; k++) rB[i * Nc + k] /= vol;
515: }
516: }
517: }
518: if (D) {
519: PetscInt i, j, k, l, m;
521: if (type == PETSC_SCALAR) {
522: PetscScalar * sD = (PetscScalar *) D;
523: PetscScalar * sqD = (PetscScalar *) qD;
525: for (i = 0; i < numPoints; i++) {
526: PetscReal vol = 0.;
528: for (j = 0; j < Nc * dimC; j++) {sD[i * Nc * dimC + j] = 0.;}
529: for (k = 0; k < Nq; k++) {
530: vol += geom->detJ[i * Nq + k] * weights[k];
531: for (j = 0; j < Nc; j++) {
532: PetscScalar pD[3] = {0.,0.,0.};
534: for (l = 0; l < dimC; l++) {
535: for (m = 0; m < dim; m++) {
536: pD[l] += geom->invJ[((i * Nq + k) * dimC + m) * dimC + l] * sqD[((i * Nq + k) * Nc + j) * dim + m];
537: }
538: }
539: for (l = 0; l < dimC; l++) {
540: sD[(i * Nc + j) * dimC + l] += geom->detJ[i * Nq + k] * weights[k] * pD[l];
541: }
542: }
543: }
544: for (k = 0; k < Nc * dimC; k++) sD[i * Nc * dimC + k] /= vol;
545: }
546: } else {
547: PetscReal * rD = (PetscReal *) D;
548: PetscReal * rqD = (PetscReal *) qD;
550: for (i = 0; i < numPoints; i++) {
551: PetscReal vol = 0.;
553: for (j = 0; j < Nc * dimC; j++) {rD[i * Nc * dimC + j] = 0.;}
554: for (k = 0; k < Nq; k++) {
555: vol += geom->detJ[i * Nq + k] * weights[k];
556: for (j = 0; j < Nc; j++) {
557: PetscReal pD[3] = {0.,0.,0.};
559: for (l = 0; l < dimC; l++) {
560: for (m = 0; m < dim; m++) {
561: pD[l] += geom->invJ[((i * Nq + k) * dimC + m) * dimC + l] * rqD[((i * Nq + k) * Nc + j) * dim + m];
562: }
563: }
564: for (l = 0; l < dimC; l++) {
565: rD[(i * Nc + j) * dimC + l] += geom->detJ[i * Nq + k] * weights[k] * pD[l];
566: }
567: }
568: }
569: for (k = 0; k < Nc * dimC; k++) rD[i * Nc * dimC + k] /= vol;
570: }
571: }
572: }
573: if (H) {
574: PetscInt i, j, k, l, m, q, r;
576: if (type == PETSC_SCALAR) {
577: PetscScalar * sH = (PetscScalar *) H;
578: PetscScalar * sqH = (PetscScalar *) qH;
580: for (i = 0; i < numPoints; i++) {
581: PetscReal vol = 0.;
583: for (j = 0; j < Nc * dimC * dimC; j++) {sH[i * Nc * dimC * dimC + j] = 0.;}
584: for (k = 0; k < Nq; k++) {
585: const PetscReal *invJ = &geom->invJ[(i * Nq + k) * dimC * dimC];
587: vol += geom->detJ[i * Nq + k] * weights[k];
588: for (j = 0; j < Nc; j++) {
589: PetscScalar pH[3][3] = {{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}};
590: const PetscScalar *spH = &sqH[((i * Nq + k) * Nc + j) * dimC * dimC];
592: for (l = 0; l < dimC; l++) {
593: for (m = 0; m < dimC; m++) {
594: for (q = 0; q < dim; q++) {
595: for (r = 0; r < dim; r++) {
596: pH[l][m] += invJ[q * dimC + l] * invJ[r * dimC + m] * spH[q * dim + r];
597: }
598: }
599: }
600: }
601: for (l = 0; l < dimC; l++) {
602: for (m = 0; m < dimC; m++) {
603: sH[(i * Nc + j) * dimC * dimC + l * dimC + m] += geom->detJ[i * Nq + k] * weights[k] * pH[l][m];
604: }
605: }
606: }
607: }
608: for (k = 0; k < Nc * dimC * dimC; k++) sH[i * Nc * dimC * dimC + k] /= vol;
609: }
610: } else {
611: PetscReal * rH = (PetscReal *) H;
612: PetscReal * rqH = (PetscReal *) qH;
614: for (i = 0; i < numPoints; i++) {
615: PetscReal vol = 0.;
617: for (j = 0; j < Nc * dimC * dimC; j++) {rH[i * Nc * dimC * dimC + j] = 0.;}
618: for (k = 0; k < Nq; k++) {
619: const PetscReal *invJ = &geom->invJ[(i * Nq + k) * dimC * dimC];
621: vol += geom->detJ[i * Nq + k] * weights[k];
622: for (j = 0; j < Nc; j++) {
623: PetscReal pH[3][3] = {{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}};
624: const PetscReal *rpH = &rqH[((i * Nq + k) * Nc + j) * dimC * dimC];
626: for (l = 0; l < dimC; l++) {
627: for (m = 0; m < dimC; m++) {
628: for (q = 0; q < dim; q++) {
629: for (r = 0; r < dim; r++) {
630: pH[l][m] += invJ[q * dimC + l] * invJ[r * dimC + m] * rpH[q * dim + r];
631: }
632: }
633: }
634: }
635: for (l = 0; l < dimC; l++) {
636: for (m = 0; m < dimC; m++) {
637: rH[(i * Nc + j) * dimC * dimC + l * dimC + m] += geom->detJ[i * Nq + k] * weights[k] * pH[l][m];
638: }
639: }
640: }
641: }
642: for (k = 0; k < Nc * dimC * dimC; k++) rH[i * Nc * dimC * dimC + k] /= vol;
643: }
644: }
645: }
646: if (B) {DMRestoreWorkArray(field->dm, N, mpitype, &qB);}
647: if (D) {DMRestoreWorkArray(field->dm, N * dimC, mpitype, &qD);}
648: if (H) {DMRestoreWorkArray(field->dm, N * dimC * dimC, mpitype, &qH);}
649: PetscFEGeomDestroy(&geom);
650: PetscQuadratureDestroy(&quad);
651: return(0);
652: }
654: static PetscErrorCode DMFieldGetDegree_DS(DMField field, IS pointIS, PetscInt *minDegree, PetscInt *maxDegree)
655: {
656: DMField_DS *dsfield;
657: PetscObject disc;
658: PetscInt h, imin, imax;
659: PetscClassId id;
663: dsfield = (DMField_DS *) field->data;
664: ISGetMinMax(pointIS,&imin,&imax);
665: if (imin >= imax) {
666: h = 0;
667: } else {
668: for (h = 0; h < dsfield->height; h++) {
669: PetscInt hEnd;
671: DMPlexGetHeightStratum(field->dm,h,NULL,&hEnd);
672: if (imin < hEnd) break;
673: }
674: }
675: DMFieldDSGetHeightDisc(field,h,&disc);
676: PetscObjectGetClassId(disc,&id);
677: if (id == PETSCFE_CLASSID) {
678: PetscFE fe = (PetscFE) disc;
679: PetscSpace sp;
681: PetscFEGetBasisSpace(fe, &sp);
682: PetscSpaceGetDegree(sp, minDegree, maxDegree);
683: }
684: return(0);
685: }
687: static PetscErrorCode DMFieldCreateDefaultQuadrature_DS(DMField field, IS pointIS, PetscQuadrature *quad)
688: {
689: PetscInt h, dim, imax, imin, cellHeight;
690: DM dm;
691: DMField_DS *dsfield;
692: PetscObject disc;
693: PetscFE fe;
694: PetscClassId id;
698: dm = field->dm;
699: dsfield = (DMField_DS *) field->data;
700: ISGetMinMax(pointIS,&imin,&imax);
701: DMGetDimension(dm,&dim);
702: for (h = 0; h <= dim; h++) {
703: PetscInt hStart, hEnd;
705: DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);
706: if (imax >= hStart && imin < hEnd) break;
707: }
708: DMPlexGetVTKCellHeight(dm, &cellHeight);
709: h -= cellHeight;
710: *quad = NULL;
711: if (h < dsfield->height) {
712: DMFieldDSGetHeightDisc(field,h,&disc);
713: PetscObjectGetClassId(disc,&id);
714: if (id != PETSCFE_CLASSID) return(0);
715: fe = (PetscFE) disc;
716: PetscFEGetQuadrature(fe,quad);
717: PetscObjectReference((PetscObject)*quad);
718: }
719: return(0);
720: }
722: static PetscErrorCode DMFieldComputeFaceData_DS(DMField field, IS pointIS, PetscQuadrature quad, PetscFEGeom *geom)
723: {
724: const PetscInt *points;
725: PetscInt p, dim, dE, numFaces, Nq;
726: PetscInt maxDegree;
727: DMLabel depthLabel;
728: IS cellIS;
729: DM dm = field->dm;
730: PetscErrorCode ierr;
733: dim = geom->dim;
734: dE = geom->dimEmbed;
735: DMPlexGetDepthLabel(dm, &depthLabel);
736: DMLabelGetStratumIS(depthLabel, dim + 1, &cellIS);
737: DMFieldGetDegree(field,cellIS,NULL,&maxDegree);
738: ISGetIndices(pointIS, &points);
739: numFaces = geom->numCells;
740: Nq = geom->numPoints;
741: /* First, set local faces and flip normals so that they are outward for the first supporting cell */
742: for (p = 0; p < numFaces; p++) {
743: PetscInt point = points[p];
744: PetscInt suppSize, s, coneSize, c, numChildren;
745: const PetscInt *supp, *cone, *ornt;
747: DMPlexGetTreeChildren(dm, point, &numChildren, NULL);
748: if (numChildren) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children");
749: DMPlexGetSupportSize(dm, point, &suppSize);
750: if (suppSize > 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D has %D support, expected at most 2\n", point, suppSize);
751: if (!suppSize) continue;
752: DMPlexGetSupport(dm, point, &supp);
753: for (s = 0; s < suppSize; ++s) {
754: DMPlexGetConeSize(dm, supp[s], &coneSize);
755: DMPlexGetCone(dm, supp[s], &cone);
756: for (c = 0; c < coneSize; ++c) if (cone[c] == point) break;
757: if (c == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid connectivity: point %D not found in cone of support point %D", point, supp[s]);
758: geom->face[p][s] = c;
759: }
760: DMPlexGetConeOrientation(dm, supp[0], &ornt);
761: if (ornt[geom->face[p][0]] < 0) {
762: PetscInt Np = geom->numPoints, q, dE = geom->dimEmbed, d;
764: for (q = 0; q < Np; ++q) for (d = 0; d < dE; ++d) geom->n[(p*Np + q)*dE + d] = -geom->n[(p*Np + q)*dE + d];
765: }
766: }
767: if (maxDegree <= 1) {
768: PetscInt numCells, offset, *cells;
769: PetscFEGeom *cellGeom;
770: IS suppIS;
772: for (p = 0, numCells = 0; p < numFaces; p++) {
773: PetscInt point = points[p];
774: PetscInt numSupp, numChildren;
776: DMPlexGetTreeChildren(dm, point, &numChildren, NULL);
777: if (numChildren) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children");
778: DMPlexGetSupportSize(dm, point,&numSupp);
779: if (numSupp > 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D has %D support, expected at most 2\n", point, numSupp);
780: numCells += numSupp;
781: }
782: PetscMalloc1(numCells, &cells);
783: for (p = 0, offset = 0; p < numFaces; p++) {
784: PetscInt point = points[p];
785: PetscInt numSupp, s;
786: const PetscInt *supp;
788: DMPlexGetSupportSize(dm, point,&numSupp);
789: DMPlexGetSupport(dm, point, &supp);
790: for (s = 0; s < numSupp; s++, offset++) {
791: cells[offset] = supp[s];
792: }
793: }
794: ISCreateGeneral(PETSC_COMM_SELF,numCells,cells,PETSC_USE_POINTER, &suppIS);
795: DMFieldCreateFEGeom(field,suppIS,quad,PETSC_FALSE,&cellGeom);
796: for (p = 0, offset = 0; p < numFaces; p++) {
797: PetscInt point = points[p];
798: PetscInt numSupp, s, q;
799: const PetscInt *supp;
801: DMPlexGetSupportSize(dm, point,&numSupp);
802: DMPlexGetSupport(dm, point, &supp);
803: for (s = 0; s < numSupp; s++, offset++) {
804: for (q = 0; q < Nq * dE * dE; q++) {
805: geom->suppJ[s][p * Nq * dE * dE + q] = cellGeom->J[offset * Nq * dE * dE + q];
806: geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q];
807: }
808: for (q = 0; q < Nq; q++) geom->suppDetJ[s][p * Nq + q] = cellGeom->detJ[offset * Nq + q];
809: }
810: }
811: PetscFEGeomDestroy(&cellGeom);
812: ISDestroy(&suppIS);
813: PetscFree(cells);
814: } else {
815: PetscObject faceDisc, cellDisc;
816: PetscClassId faceId, cellId;
817: PetscDualSpace dsp;
818: DM K;
819: DMPolytopeType ct;
820: PetscInt (*co)[2][3];
821: PetscInt coneSize;
822: PetscInt **counts;
823: PetscInt f, i, o, q, s;
824: const PetscInt *coneK;
825: PetscInt eStart, minOrient, maxOrient, numOrient;
826: PetscInt *orients;
827: PetscReal **orientPoints;
828: PetscReal *cellPoints;
829: PetscReal *dummyWeights;
830: PetscQuadrature cellQuad = NULL;
832: DMFieldDSGetHeightDisc(field, 1, &faceDisc);
833: DMFieldDSGetHeightDisc(field, 0, &cellDisc);
834: PetscObjectGetClassId(faceDisc,&faceId);
835: PetscObjectGetClassId(cellDisc,&cellId);
836: if (faceId != PETSCFE_CLASSID || cellId != PETSCFE_CLASSID) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Not supported\n");
837: PetscFEGetDualSpace((PetscFE)cellDisc, &dsp);
838: PetscDualSpaceGetDM(dsp, &K);
839: DMPlexGetHeightStratum(K, 1, &eStart, NULL);
840: DMPlexGetCellType(K, eStart, &ct);
841: DMPlexGetConeSize(K,0,&coneSize);
842: DMPlexGetCone(K,0,&coneK);
843: PetscMalloc2(numFaces, &co, coneSize, &counts);
844: PetscMalloc1(dE*Nq, &cellPoints);
845: PetscMalloc1(Nq, &dummyWeights);
846: PetscQuadratureCreate(PETSC_COMM_SELF, &cellQuad);
847: PetscQuadratureSetData(cellQuad, dE, 1, Nq, cellPoints, dummyWeights);
848: minOrient = PETSC_MAX_INT;
849: maxOrient = PETSC_MIN_INT;
850: for (p = 0; p < numFaces; p++) { /* record the orientation of the facet wrt the support cells */
851: PetscInt point = points[p];
852: PetscInt numSupp, numChildren;
853: const PetscInt *supp;
855: DMPlexGetTreeChildren(dm, point, &numChildren, NULL);
856: if (numChildren) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face data not valid for facets with children");
857: DMPlexGetSupportSize(dm, point,&numSupp);
858: if (numSupp > 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D has %D support, expected at most 2\n", point, numSupp);
859: DMPlexGetSupport(dm, point, &supp);
860: for (s = 0; s < numSupp; s++) {
861: PetscInt cell = supp[s];
862: PetscInt numCone;
863: const PetscInt *cone, *orient;
865: DMPlexGetConeSize(dm, cell, &numCone);
866: if (numCone != coneSize) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Support point does not match reference element");
867: DMPlexGetCone(dm, cell, &cone);
868: DMPlexGetConeOrientation(dm, cell, &orient);
869: for (f = 0; f < coneSize; f++) {
870: if (cone[f] == point) break;
871: }
872: co[p][s][0] = f;
873: co[p][s][1] = orient[f];
874: co[p][s][2] = cell;
875: minOrient = PetscMin(minOrient, orient[f]);
876: maxOrient = PetscMax(maxOrient, orient[f]);
877: }
878: for (; s < 2; s++) {
879: co[p][s][0] = -1;
880: co[p][s][1] = -1;
881: co[p][s][2] = -1;
882: }
883: }
884: numOrient = maxOrient + 1 - minOrient;
885: DMPlexGetCone(K,0,&coneK);
886: /* count all (face,orientation) doubles that appear */
887: PetscCalloc2(numOrient,&orients,numOrient,&orientPoints);
888: for (f = 0; f < coneSize; f++) {PetscCalloc1(numOrient+1, &counts[f]);}
889: for (p = 0; p < numFaces; p++) {
890: for (s = 0; s < 2; s++) {
891: if (co[p][s][0] >= 0) {
892: counts[co[p][s][0]][co[p][s][1] - minOrient]++;
893: orients[co[p][s][1] - minOrient]++;
894: }
895: }
896: }
897: for (o = 0; o < numOrient; o++) {
898: if (orients[o]) {
899: PetscInt orient = o + minOrient;
900: PetscInt q;
902: PetscMalloc1(Nq * dim, &orientPoints[o]);
903: /* rotate the quadrature points appropriately */
904: switch (ct) {
905: case DM_POLYTOPE_POINT: break;
906: case DM_POLYTOPE_SEGMENT:
907: if (orient == -2 || orient == 1) {
908: for (q = 0; q < Nq; q++) {
909: orientPoints[o][q] = -geom->xi[q];
910: }
911: } else {
912: for (q = 0; q < Nq; q++) {
913: orientPoints[o][q] = geom->xi[q];
914: }
915: }
916: break;
917: case DM_POLYTOPE_TRIANGLE:
918: for (q = 0; q < Nq; q++) {
919: PetscReal lambda[3];
920: PetscReal lambdao[3];
922: /* convert to barycentric */
923: lambda[0] = - (geom->xi[2 * q] + geom->xi[2 * q + 1]) / 2.;
924: lambda[1] = (geom->xi[2 * q] + 1.) / 2.;
925: lambda[2] = (geom->xi[2 * q + 1] + 1.) / 2.;
926: if (orient >= 0) {
927: for (i = 0; i < 3; i++) {
928: lambdao[i] = lambda[(orient + i) % 3];
929: }
930: } else {
931: for (i = 0; i < 3; i++) {
932: lambdao[i] = lambda[(-(orient + i) + 3) % 3];
933: }
934: }
935: /* convert to coordinates */
936: orientPoints[o][2 * q + 0] = -(lambdao[0] + lambdao[2]) + lambdao[1];
937: orientPoints[o][2 * q + 1] = -(lambdao[0] + lambdao[1]) + lambdao[2];
938: }
939: break;
940: case DM_POLYTOPE_QUADRILATERAL:
941: for (q = 0; q < Nq; q++) {
942: PetscReal xi[2], xio[2];
943: PetscInt oabs = (orient >= 0) ? orient : -(orient + 1);
945: xi[0] = geom->xi[2 * q];
946: xi[1] = geom->xi[2 * q + 1];
947: switch (oabs) {
948: case 1:
949: xio[0] = xi[1];
950: xio[1] = -xi[0];
951: break;
952: case 2:
953: xio[0] = -xi[0];
954: xio[1] = -xi[1];
955: case 3:
956: xio[0] = -xi[1];
957: xio[1] = xi[0];
958: case 0:
959: default:
960: xio[0] = xi[0];
961: xio[1] = xi[1];
962: break;
963: }
964: if (orient < 0) {
965: xio[0] = -xio[0];
966: }
967: orientPoints[o][2 * q + 0] = xio[0];
968: orientPoints[o][2 * q + 1] = xio[1];
969: }
970: break;
971: default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cell type %s not yet supported\n", DMPolytopeTypes[ct]);
972: }
973: }
974: }
975: for (f = 0; f < coneSize; f++) {
976: PetscInt face = coneK[f];
977: PetscReal v0[3];
978: PetscReal J[9], detJ;
979: PetscInt numCells, offset;
980: PetscInt *cells;
981: IS suppIS;
983: DMPlexComputeCellGeometryFEM(K, face, NULL, v0, J, NULL, &detJ);
984: for (o = 0; o <= numOrient; o++) {
985: PetscFEGeom *cellGeom;
987: if (!counts[f][o]) continue;
988: /* If this (face,orientation) double appears,
989: * convert the face quadrature points into volume quadrature points */
990: for (q = 0; q < Nq; q++) {
991: PetscReal xi0[3] = {-1., -1., -1.};
993: CoordinatesRefToReal(dE, dim, xi0, v0, J, &orientPoints[o][dim * q + 0], &cellPoints[dE * q + 0]);
994: }
995: for (p = 0, numCells = 0; p < numFaces; p++) {
996: for (s = 0; s < 2; s++) {
997: if (co[p][s][0] == f && co[p][s][1] == o + minOrient) numCells++;
998: }
999: }
1000: PetscMalloc1(numCells, &cells);
1001: for (p = 0, offset = 0; p < numFaces; p++) {
1002: for (s = 0; s < 2; s++) {
1003: if (co[p][s][0] == f && co[p][s][1] == o + minOrient) {
1004: cells[offset++] = co[p][s][2];
1005: }
1006: }
1007: }
1008: ISCreateGeneral(PETSC_COMM_SELF,numCells,cells,PETSC_USE_POINTER, &suppIS);
1009: DMFieldCreateFEGeom(field,suppIS,cellQuad,PETSC_FALSE,&cellGeom);
1010: for (p = 0, offset = 0; p < numFaces; p++) {
1011: for (s = 0; s < 2; s++) {
1012: if (co[p][s][0] == f && co[p][s][1] == o + minOrient) {
1013: for (q = 0; q < Nq * dE * dE; q++) {
1014: geom->suppJ[s][p * Nq * dE * dE + q] = cellGeom->J[offset * Nq * dE * dE + q];
1015: geom->suppInvJ[s][p * Nq * dE * dE + q] = cellGeom->invJ[offset * Nq * dE * dE + q];
1016: }
1017: for (q = 0; q < Nq; q++) geom->suppDetJ[s][p * Nq + q] = cellGeom->detJ[offset * Nq + q];
1018: offset++;
1019: }
1020: }
1021: }
1022: PetscFEGeomDestroy(&cellGeom);
1023: ISDestroy(&suppIS);
1024: PetscFree(cells);
1025: }
1026: }
1027: for (o = 0; o < numOrient; o++) {
1028: if (orients[o]) {
1029: PetscFree(orientPoints[o]);
1030: }
1031: }
1032: PetscFree2(orients,orientPoints);
1033: PetscQuadratureDestroy(&cellQuad);
1034: for (f = 0; f < coneSize; f++) {PetscFree(counts[f]);}
1035: PetscFree2(co,counts);
1036: }
1037: ISRestoreIndices(pointIS, &points);
1038: ISDestroy(&cellIS);
1039: return(0);
1040: }
1042: static PetscErrorCode DMFieldInitialize_DS(DMField field)
1043: {
1045: field->ops->destroy = DMFieldDestroy_DS;
1046: field->ops->evaluate = DMFieldEvaluate_DS;
1047: field->ops->evaluateFE = DMFieldEvaluateFE_DS;
1048: field->ops->evaluateFV = DMFieldEvaluateFV_DS;
1049: field->ops->getDegree = DMFieldGetDegree_DS;
1050: field->ops->createDefaultQuadrature = DMFieldCreateDefaultQuadrature_DS;
1051: field->ops->view = DMFieldView_DS;
1052: field->ops->computeFaceData = DMFieldComputeFaceData_DS;
1053: return(0);
1054: }
1056: PETSC_INTERN PetscErrorCode DMFieldCreate_DS(DMField field)
1057: {
1058: DMField_DS *dsfield;
1062: PetscNewLog(field,&dsfield);
1063: field->data = dsfield;
1064: DMFieldInitialize_DS(field);
1065: return(0);
1066: }
1068: PetscErrorCode DMFieldCreateDS(DM dm, PetscInt fieldNum, Vec vec,DMField *field)
1069: {
1070: DMField b;
1071: DMField_DS *dsfield;
1072: PetscObject disc = NULL;
1073: PetscBool isContainer = PETSC_FALSE;
1074: PetscClassId id = -1;
1075: PetscInt numComponents = -1, dsNumFields;
1076: PetscSection section;
1080: DMGetLocalSection(dm,§ion);
1081: PetscSectionGetFieldComponents(section,fieldNum,&numComponents);
1082: DMGetNumFields(dm,&dsNumFields);
1083: if (dsNumFields) {DMGetField(dm,fieldNum,NULL,&disc);}
1084: if (disc) {
1085: PetscObjectGetClassId(disc,&id);
1086: isContainer = (id == PETSC_CONTAINER_CLASSID) ? PETSC_TRUE : PETSC_FALSE;
1087: }
1088: if (!disc || isContainer) {
1089: MPI_Comm comm = PetscObjectComm((PetscObject) dm);
1090: PetscInt cStart, cEnd, dim, cellHeight;
1091: PetscInt localConeSize = 0, coneSize;
1092: PetscFE fe;
1093: PetscDualSpace Q;
1094: PetscSpace P;
1095: DM K;
1096: PetscQuadrature quad, fquad;
1097: PetscBool isSimplex;
1099: DMPlexGetVTKCellHeight(dm, &cellHeight);
1100: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
1101: DMGetDimension(dm, &dim);
1102: if (cEnd > cStart) {
1103: DMPlexGetConeSize(dm, cStart, &localConeSize);
1104: }
1105: MPI_Allreduce(&localConeSize,&coneSize,1,MPIU_INT,MPI_MAX,comm);
1106: isSimplex = (coneSize == (dim + 1)) ? PETSC_TRUE : PETSC_FALSE;
1107: PetscSpaceCreate(PETSC_COMM_SELF, &P);
1108: PetscSpaceSetType(P,PETSCSPACEPOLYNOMIAL);
1109: PetscSpaceSetDegree(P, 1, PETSC_DETERMINE);
1110: PetscSpaceSetNumComponents(P, numComponents);
1111: PetscSpaceSetNumVariables(P, dim);
1112: PetscSpacePolynomialSetTensor(P, isSimplex ? PETSC_FALSE : PETSC_TRUE);
1113: PetscSpaceSetUp(P);
1114: PetscDualSpaceCreate(PETSC_COMM_SELF, &Q);
1115: PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);
1116: PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1117: PetscDualSpaceSetDM(Q, K);
1118: DMDestroy(&K);
1119: PetscDualSpaceSetNumComponents(Q, numComponents);
1120: PetscDualSpaceSetOrder(Q, 1);
1121: PetscDualSpaceLagrangeSetTensor(Q, isSimplex ? PETSC_FALSE : PETSC_TRUE);
1122: PetscDualSpaceSetUp(Q);
1123: PetscFECreate(PETSC_COMM_SELF, &fe);
1124: PetscFESetType(fe,PETSCFEBASIC);
1125: PetscFESetBasisSpace(fe, P);
1126: PetscFESetDualSpace(fe, Q);
1127: PetscFESetNumComponents(fe, numComponents);
1128: PetscFESetUp(fe);
1129: PetscSpaceDestroy(&P);
1130: PetscDualSpaceDestroy(&Q);
1131: if (isSimplex) {
1132: PetscDTStroudConicalQuadrature(dim, 1, 1, -1.0, 1.0, &quad);
1133: PetscDTStroudConicalQuadrature(dim-1, 1, 1, -1.0, 1.0, &fquad);
1134: }
1135: else {
1136: PetscDTGaussTensorQuadrature(dim, 1, 1, -1.0, 1.0, &quad);
1137: PetscDTGaussTensorQuadrature(dim-1, 1, 1, -1.0, 1.0, &fquad);
1138: }
1139: PetscFESetQuadrature(fe, quad);
1140: PetscFESetFaceQuadrature(fe, fquad);
1141: PetscQuadratureDestroy(&quad);
1142: PetscQuadratureDestroy(&fquad);
1143: disc = (PetscObject) fe;
1144: } else {
1145: PetscObjectReference(disc);
1146: }
1147: PetscObjectGetClassId(disc,&id);
1148: if (id == PETSCFE_CLASSID) {
1149: PetscFE fe = (PetscFE) disc;
1151: PetscFEGetNumComponents(fe,&numComponents);
1152: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented");
1153: DMFieldCreate(dm,numComponents,DMFIELD_VERTEX,&b);
1154: DMFieldSetType(b,DMFIELDDS);
1155: dsfield = (DMField_DS *) b->data;
1156: dsfield->fieldNum = fieldNum;
1157: DMGetDimension(dm,&dsfield->height);
1158: dsfield->height++;
1159: PetscCalloc1(dsfield->height,&dsfield->disc);
1160: dsfield->disc[0] = disc;
1161: PetscObjectReference((PetscObject)vec);
1162: dsfield->vec = vec;
1163: *field = b;
1165: return(0);
1166: }