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