Actual source code: dmfield.c
1: #include <petsc/private/dmfieldimpl.h>
2: #include <petsc/private/petscfeimpl.h>
3: #include <petscdmplex.h>
5: const char *const DMFieldContinuities[] = {
6: "VERTEX",
7: "EDGE",
8: "FACET",
9: "CELL",
10: NULL
11: };
13: PETSC_INTERN PetscErrorCode DMFieldCreate(DM dm,PetscInt numComponents,DMFieldContinuity continuity,DMField *field)
14: {
16: DMField b;
21: DMFieldInitializePackage();
23: PetscHeaderCreate(b,DMFIELD_CLASSID,"DMField","Field over DM","DM",PetscObjectComm((PetscObject)dm),DMFieldDestroy,DMFieldView);
24: PetscObjectReference((PetscObject)dm);
25: b->dm = dm;
26: b->continuity = continuity;
27: b->numComponents = numComponents;
28: *field = b;
29: return(0);
30: }
32: /*@
33: DMFieldDestroy - destroy a DMField
35: Collective
37: Input Parameter:
38: . field - address of DMField
40: Level: advanced
42: .seealso: DMFieldCreate()
43: @*/
44: PetscErrorCode DMFieldDestroy(DMField *field)
45: {
49: if (!*field) return(0);
51: if (--((PetscObject)(*field))->refct > 0) {*field = NULL; return(0);}
52: if ((*field)->ops->destroy) {(*(*field)->ops->destroy)(*field);}
53: DMDestroy(&((*field)->dm));
54: PetscHeaderDestroy(field);
55: return(0);
56: }
58: /*@C
59: DMFieldView - view a DMField
61: Collective
63: Input Parameters:
64: + field - DMField
65: - viewer - viewer to display field, for example PETSC_VIEWER_STDOUT_WORLD
67: Level: advanced
69: .seealso: DMFieldCreate()
70: @*/
71: PetscErrorCode DMFieldView(DMField field,PetscViewer viewer)
72: {
73: PetscErrorCode ierr;
74: PetscBool iascii;
78: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)field),&viewer);}
81: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
82: if (iascii) {
83: PetscObjectPrintClassNamePrefixType((PetscObject)field,viewer);
84: PetscViewerASCIIPushTab(viewer);
85: PetscViewerASCIIPrintf(viewer,"%D components\n",field->numComponents);
86: PetscViewerASCIIPrintf(viewer,"%s continuity\n",DMFieldContinuities[field->continuity]);
87: PetscViewerPushFormat(viewer,PETSC_VIEWER_DEFAULT);
88: DMView(field->dm,viewer);
89: PetscViewerPopFormat(viewer);
90: }
91: if (field->ops->view) {(*field->ops->view)(field,viewer);}
92: if (iascii) {
93: PetscViewerASCIIPopTab(viewer);
94: }
95: return(0);
96: }
98: /*@C
99: DMFieldSetType - set the DMField implementation
101: Collective on field
103: Input Parameters:
104: + field - the DMField context
105: - type - a known method
107: Notes:
108: See "include/petscvec.h" for available methods (for instance)
109: + DMFIELDDA - a field defined only by its values at the corners of a DMDA
110: . DMFIELDDS - a field defined by a discretization over a mesh set with DMSetField()
111: - DMFIELDSHELL - a field defined by arbitrary callbacks
113: Level: advanced
115: .seealso: DMFieldType,
116: @*/
117: PetscErrorCode DMFieldSetType(DMField field,DMFieldType type)
118: {
119: PetscErrorCode ierr,(*r)(DMField);
120: PetscBool match;
126: PetscObjectTypeCompare((PetscObject)field,type,&match);
127: if (match) return(0);
129: PetscFunctionListFind(DMFieldList,type,&r);
130: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested DMField type %s",type);
131: /* Destroy the previous private DMField context */
132: if (field->ops->destroy) {
133: (*(field)->ops->destroy)(field);
134: }
135: PetscMemzero(field->ops,sizeof(*field->ops));
136: PetscObjectChangeTypeName((PetscObject)field,type);
137: field->ops->create = r;
138: (*r)(field);
139: return(0);
140: }
142: /*@C
143: DMFieldGetType - Gets the DMField type name (as a string) from the DMField.
145: Not Collective
147: Input Parameter:
148: . field - The DMField context
150: Output Parameter:
151: . type - The DMField type name
153: Level: advanced
155: .seealso: DMFieldSetType()
156: @*/
157: PetscErrorCode DMFieldGetType(DMField field, DMFieldType *type)
158: {
164: DMFieldRegisterAll();
165: *type = ((PetscObject)field)->type_name;
166: return(0);
167: }
169: /*@
170: DMFieldGetNumComponents - Returns the number of components in the field
172: Not collective
174: Input Parameter:
175: . field - The DMField object
177: Output Parameter:
178: . nc - The number of field components
180: Level: intermediate
182: .seealso: DMFieldEvaluate()
183: @*/
184: PetscErrorCode DMFieldGetNumComponents(DMField field, PetscInt *nc)
185: {
189: *nc = field->numComponents;
190: return(0);
191: }
193: /*@
194: DMFieldGetDM - Returns the DM for the manifold over which the field is defined.
196: Not collective
198: Input Parameter:
199: . field - The DMField object
201: Output Parameter:
202: . dm - The DM object
204: Level: intermediate
206: .seealso: DMFieldEvaluate()
207: @*/
208: PetscErrorCode DMFieldGetDM(DMField field, DM *dm)
209: {
213: *dm = field->dm;
214: return(0);
215: }
217: /*@
218: DMFieldEvaluate - Evaluate the field and its derivatives on a set of points
220: Collective on points
222: Input Parameters:
223: + field - The DMField object
224: . points - The points at which to evaluate the field. Should have size d x n,
225: where d is the coordinate dimension of the manifold and n is the number
226: of points
227: - datatype - The PetscDataType of the output arrays: either PETSC_REAL or PETSC_SCALAR.
228: If the field is complex and datatype is PETSC_REAL, the real part of the
229: field is returned.
231: Output Parameters:
232: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
233: If B is not NULL, the values of the field are written in this array, varying first by component,
234: then by point.
235: . D - pointer to data of size d * c * n * sizeof(datatype).
236: If D is not NULL, the values of the field's spatial derivatives are written in this array,
237: varying first by the partial derivative component, then by field component, then by point.
238: - H - pointer to data of size d * d * c * n * sizeof(datatype).
239: If H is not NULL, the values of the field's second spatial derivatives are written in this array,
240: varying first by the second partial derivative component, then by field component, then by point.
242: Level: intermediate
244: .seealso: DMFieldGetDM(), DMFieldGetNumComponents(), DMFieldEvaluateFE(), DMFieldEvaluateFV()
245: @*/
246: PetscErrorCode DMFieldEvaluate(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
247: {
256: if (field->ops->evaluate) {
257: (*field->ops->evaluate) (field, points, datatype, B, D, H);
258: } else SETERRQ (PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented for this type");
259: return(0);
260: }
262: /*@
263: DMFieldEvaluateFE - Evaluate the field and its derivatives on a set of points mapped from
264: quadrature points on a reference point. The derivatives are taken with respect to the
265: reference coordinates.
267: Not collective
269: Input Parameters:
270: + field - The DMField object
271: . cellIS - Index set for cells on which to evaluate the field
272: . points - The quadature containing the points in the reference cell at which to evaluate the field.
273: - datatype - The PetscDataType of the output arrays: either PETSC_REAL or PETSC_SCALAR.
274: If the field is complex and datatype is PETSC_REAL, the real part of the
275: field is returned.
277: Output Parameters:
278: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
279: If B is not NULL, the values of the field are written in this array, varying first by component,
280: then by point.
281: . D - pointer to data of size d * c * n * sizeof(datatype).
282: If D is not NULL, the values of the field's spatial derivatives are written in this array,
283: varying first by the partial derivative component, then by field component, then by point.
284: - H - pointer to data of size d * d * c * n * sizeof(datatype).
285: If H is not NULL, the values of the field's second spatial derivatives are written in this array,
286: varying first by the second partial derivative component, then by field component, then by point.
288: Level: intermediate
290: .seealso: DMFieldGetNumComponents(), DMFieldEvaluate(), DMFieldEvaluateFV()
291: @*/
292: PetscErrorCode DMFieldEvaluateFE(DMField field, IS cellIS, PetscQuadrature points, PetscDataType datatype, void *B, void *D, void *H)
293: {
303: if (field->ops->evaluateFE) {
304: (*field->ops->evaluateFE) (field, cellIS, points, datatype, B, D, H);
305: } else SETERRQ (PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented for this type");
306: return(0);
307: }
309: /*@
310: DMFieldEvaluateFV - Evaluate the mean of a field and its finite volume derivatives on a set of points.
312: Not collective
314: Input Parameters:
315: + field - The DMField object
316: . cellIS - Index set for cells on which to evaluate the field
317: - datatype - The PetscDataType of the output arrays: either PETSC_REAL or PETSC_SCALAR.
318: If the field is complex and datatype is PETSC_REAL, the real part of the
319: field is returned.
321: Output Parameters:
322: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
323: If B is not NULL, the values of the field are written in this array, varying first by component,
324: then by point.
325: . D - pointer to data of size d * c * n * sizeof(datatype).
326: If D is not NULL, the values of the field's spatial derivatives are written in this array,
327: varying first by the partial derivative component, then by field component, then by point.
328: - H - pointer to data of size d * d * c * n * sizeof(datatype).
329: If H is not NULL, the values of the field's second spatial derivatives are written in this array,
330: varying first by the second partial derivative component, then by field component, then by point.
332: Level: intermediate
334: .seealso: DMFieldGetNumComponents(), DMFieldEvaluate(), DMFieldEvaluateFE()
335: @*/
336: PetscErrorCode DMFieldEvaluateFV(DMField field, IS cellIS, PetscDataType datatype, void *B, void *D, void *H)
337: {
346: if (field->ops->evaluateFV) {
347: (*field->ops->evaluateFV) (field, cellIS, datatype, B, D, H);
348: } else SETERRQ (PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented for this type");
349: return(0);
350: }
352: /*@
353: DMFieldGetDegree - Get the polynomial degree of a field when pulled back onto the
354: reference element
356: Not collective
358: Input Parameters:
359: + field - the DMField object
360: - cellIS - the index set of points over which we want know the invariance
362: Output Parameters:
363: + minDegree - the degree of the largest polynomial space contained in the field on each element
364: - maxDegree - the largest degree of the smallest polynomial space containing the field on any element
366: Level: intermediate
368: .seealso: DMFieldEvaluateFE()
369: @*/
370: PetscErrorCode DMFieldGetDegree(DMField field, IS cellIS, PetscInt *minDegree, PetscInt *maxDegree)
371: {
380: if (minDegree) *minDegree = -1;
381: if (maxDegree) *maxDegree = PETSC_MAX_INT;
383: if (field->ops->getDegree) {
384: (*field->ops->getDegree) (field,cellIS,minDegree,maxDegree);
385: }
386: return(0);
387: }
389: /*@
390: DMFieldCreateDefaultQuadrature - Creates a quadrature sufficient to integrate the field on the selected
391: points via pullback onto the reference element
393: Not collective
395: Input Parameters:
396: + field - the DMField object
397: - pointIS - the index set of points over which we wish to integrate the field
399: Output Parameter:
400: . quad - a PetscQuadrature object
402: Level: developer
404: .seealso: DMFieldEvaluteFE(), DMFieldGetDegree()
405: @*/
406: PetscErrorCode DMFieldCreateDefaultQuadrature(DMField field, IS pointIS, PetscQuadrature *quad)
407: {
415: *quad = NULL;
416: if (field->ops->createDefaultQuadrature) {
417: (*field->ops->createDefaultQuadrature)(field, pointIS, quad);
418: }
419: return(0);
420: }
422: /*@C
423: DMFieldCreateFEGeom - Compute and create the geometric factors of a coordinate field
425: Not collective
427: Input Parameters:
428: + field - the DMField object
429: . pointIS - the index set of points over which we wish to integrate the field
430: . quad - the quadrature points at which to evaluate the geometric factors
431: - faceData - whether additional data for facets (the normal vectors and adjacent cells) should
432: be calculated
434: Output Parameter:
435: . geom - the geometric factors
437: Level: developer
439: .seealso: DMFieldEvaluateFE(), DMFieldCreateDefaulteQuadrature(), DMFieldGetDegree()
440: @*/
441: PetscErrorCode DMFieldCreateFEGeom(DMField field, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
442: {
443: PetscInt dim, dE;
444: PetscInt nPoints;
445: PetscInt maxDegree;
446: PetscFEGeom *g;
453: ISGetLocalSize(pointIS,&nPoints);
454: dE = field->numComponents;
455: PetscFEGeomCreate(quad,nPoints,dE,faceData,&g);
456: DMFieldEvaluateFE(field,pointIS,quad,PETSC_REAL,g->v,g->J,NULL);
457: dim = g->dim;
458: if (dE > dim) {
459: /* space out J and make square Jacobians */
460: PetscInt i, j, k, N = g->numPoints * g->numCells;
462: for (i = N-1; i >= 0; i--) {
463: PetscReal J[9];
465: for (j = 0; j < dE; j++) {
466: for (k = 0; k < dim; k++) {
467: J[j*dE + k] = g->J[i*dE*dim + j*dim + k];
468: }
469: }
470: switch (dim) {
471: case 0:
472: for (j = 0; j < dE; j++) {
473: for (k = 0; k < dE; k++) {
474: J[j * dE + k] = (j == k) ? 1. : 0.;
475: }
476: }
477: break;
478: case 1:
479: if (dE == 2) {
480: PetscReal norm = PetscSqrtReal(J[0] * J[0] + J[2] * J[2]);
482: J[1] = -J[2] / norm;
483: J[3] = J[0] / norm;
484: } else {
485: PetscReal inorm = 1./PetscSqrtReal(J[0] * J[0] + J[3] * J[3] + J[6] * J[6]);
486: PetscReal x = J[0] * inorm;
487: PetscReal y = J[3] * inorm;
488: PetscReal z = J[6] * inorm;
490: if (x > 0.) {
491: PetscReal inv1pX = 1./ (1. + x);
493: J[1] = -y; J[2] = -z;
494: J[4] = 1. - y*y*inv1pX; J[5] = -y*z*inv1pX;
495: J[7] = -y*z*inv1pX; J[8] = 1. - z*z*inv1pX;
496: } else {
497: PetscReal inv1mX = 1./ (1. - x);
499: J[1] = z; J[2] = y;
500: J[4] = -y*z*inv1mX; J[5] = 1. - y*y*inv1mX;
501: J[7] = 1. - z*z*inv1mX; J[8] = -y*z*inv1mX;
502: }
503: }
504: break;
505: case 2:
506: {
507: PetscReal inorm;
509: J[2] = J[3] * J[7] - J[6] * J[4];
510: J[5] = J[6] * J[1] - J[0] * J[7];
511: J[8] = J[0] * J[4] - J[3] * J[1];
513: inorm = 1./ PetscSqrtReal(J[2]*J[2] + J[5]*J[5] + J[8]*J[8]);
515: J[2] *= inorm;
516: J[5] *= inorm;
517: J[8] *= inorm;
518: }
519: break;
520: }
521: for (j = 0; j < dE*dE; j++) {
522: g->J[i*dE*dE + j] = J[j];
523: }
524: }
525: }
526: PetscFEGeomComplete(g);
527: DMFieldGetDegree(field,pointIS,NULL,&maxDegree);
528: g->isAffine = (maxDegree <= 1) ? PETSC_TRUE : PETSC_FALSE;
529: if (faceData) {
530: if (!field->ops->computeFaceData) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "DMField implementation does not compute face data\n");
531: (*field->ops->computeFaceData) (field, pointIS, quad, g);
532: }
533: *geom = g;
534: return(0);
535: }