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