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