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[] = {"VERTEX", "EDGE", "FACET", "CELL", NULL};
7: PETSC_INTERN PetscErrorCode DMFieldCreate(DM dm, PetscInt numComponents, DMFieldContinuity continuity, DMField *field)
8: {
9: DMField b;
13: DMFieldInitializePackage();
15: PetscHeaderCreate(b, DMFIELD_CLASSID, "DMField", "Field over DM", "DM", PetscObjectComm((PetscObject)dm), DMFieldDestroy, DMFieldView);
16: PetscObjectReference((PetscObject)dm);
17: b->dm = dm;
18: b->continuity = continuity;
19: b->numComponents = numComponents;
20: *field = b;
21: return 0;
22: }
24: /*@
25: DMFieldDestroy - destroy a `DMField`
27: Collective
29: Input Parameter:
30: . field - address of `DMField`
32: Level: advanced
34: .seealso: `DMField`, `DMFieldCreate()`
35: @*/
36: PetscErrorCode DMFieldDestroy(DMField *field)
37: {
38: if (!*field) return 0;
40: if (--((PetscObject)(*field))->refct > 0) {
41: *field = NULL;
42: return 0;
43: }
44: PetscTryTypeMethod((*field), destroy);
45: DMDestroy(&((*field)->dm));
46: PetscHeaderDestroy(field);
47: return 0;
48: }
50: /*@C
51: DMFieldView - view a DMField
53: Collective
55: Input Parameters:
56: + field - DMField
57: - viewer - viewer to display field, for example PETSC_VIEWER_STDOUT_WORLD
59: Level: advanced
61: .seealso: `DMField`, `DMFieldCreate()`
62: @*/
63: PetscErrorCode DMFieldView(DMField field, PetscViewer viewer)
64: {
65: PetscBool iascii;
68: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)field), &viewer);
71: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
72: if (iascii) {
73: PetscObjectPrintClassNamePrefixType((PetscObject)field, viewer);
74: PetscViewerASCIIPushTab(viewer);
75: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " components\n", field->numComponents);
76: PetscViewerASCIIPrintf(viewer, "%s continuity\n", DMFieldContinuities[field->continuity]);
77: PetscViewerPushFormat(viewer, PETSC_VIEWER_DEFAULT);
78: DMView(field->dm, viewer);
79: PetscViewerPopFormat(viewer);
80: }
81: PetscTryTypeMethod(field, view, viewer);
82: if (iascii) PetscViewerASCIIPopTab(viewer);
83: return 0;
84: }
86: /*@C
87: DMFieldSetType - set the `DMField` implementation
89: Collective on field
91: Input Parameters:
92: + field - the `DMField` context
93: - type - a known method
95: Notes:
96: See "include/petscvec.h" for available methods (for instance)
97: + `DMFIELDDA` - a field defined only by its values at the corners of a `DMDA`
98: . `DMFIELDDS` - a field defined by a discretization over a mesh set with `DMSetField()`
99: - `DMFIELDSHELL` - a field defined by arbitrary callbacks
101: Level: advanced
103: .seealso: `DMField`, `DMFieldGetType()`, `DMFieldType`,
104: @*/
105: PetscErrorCode DMFieldSetType(DMField field, DMFieldType type)
106: {
107: PetscBool match;
108: PetscErrorCode (*r)(DMField);
113: PetscObjectTypeCompare((PetscObject)field, type, &match);
114: if (match) return 0;
116: PetscFunctionListFind(DMFieldList, type, &r);
118: /* Destroy the previous private DMField context */
119: PetscTryTypeMethod(field, destroy);
121: PetscMemzero(field->ops, sizeof(*field->ops));
122: PetscObjectChangeTypeName((PetscObject)field, type);
123: field->ops->create = r;
124: (*r)(field);
125: return 0;
126: }
128: /*@C
129: DMFieldGetType - Gets the `DMFieldType` name (as a string) from the `DMField`.
131: Not Collective
133: Input Parameter:
134: . field - The `DMField` context
136: Output Parameter:
137: . type - The `DMFieldType` name
139: Level: advanced
141: .seealso: `DMField`, `DMFieldSetType()`, `DMFieldType`
142: @*/
143: PetscErrorCode DMFieldGetType(DMField field, DMFieldType *type)
144: {
147: DMFieldRegisterAll();
148: *type = ((PetscObject)field)->type_name;
149: return 0;
150: }
152: /*@
153: DMFieldGetNumComponents - Returns the number of components in the field
155: Not collective
157: Input Parameter:
158: . field - The `DMField` object
160: Output Parameter:
161: . nc - The number of field components
163: Level: intermediate
165: .seealso: `DMField`, `DMFieldEvaluate()`
166: @*/
167: PetscErrorCode DMFieldGetNumComponents(DMField field, PetscInt *nc)
168: {
171: *nc = field->numComponents;
172: return 0;
173: }
175: /*@
176: DMFieldGetDM - Returns the `DM` for the manifold over which the field is defined.
178: Not collective
180: Input Parameter:
181: . field - The `DMField` object
183: Output Parameter:
184: . dm - The `DM` object
186: Level: intermediate
188: .seealso: `DMField`, `DM`, `DMFieldEvaluate()`
189: @*/
190: PetscErrorCode DMFieldGetDM(DMField field, DM *dm)
191: {
194: *dm = field->dm;
195: return 0;
196: }
198: /*@
199: DMFieldEvaluate - Evaluate the field and its derivatives on a set of points
201: Collective on points
203: Input Parameters:
204: + field - The `DMField` object
205: . points - The points at which to evaluate the field. Should have size d x n,
206: where d is the coordinate dimension of the manifold and n is the number
207: of points
208: - datatype - The PetscDataType of the output arrays: either `PETSC_REAL` or `PETSC_SCALAR`.
209: If the field is complex and datatype is `PETSC_REAL`, the real part of the
210: field is returned.
212: Output Parameters:
213: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
214: If B is not NULL, the values of the field are written in this array, varying first by component,
215: then by point.
216: . D - pointer to data of size d * c * n * sizeof(datatype).
217: If D is not NULL, the values of the field's spatial derivatives are written in this array,
218: varying first by the partial derivative component, then by field component, then by point.
219: - H - pointer to data of size d * d * c * n * sizeof(datatype).
220: If H is not NULL, the values of the field's second spatial derivatives are written in this array,
221: varying first by the second partial derivative component, then by field component, then by point.
223: Level: intermediate
225: .seealso: `DMField`, `DMFieldGetDM()`, `DMFieldGetNumComponents()`, `DMFieldEvaluateFE()`, `DMFieldEvaluateFV()`, `PetscDataType`
226: @*/
227: PetscErrorCode DMFieldEvaluate(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
228: {
234: if (field->ops->evaluate) {
235: (*field->ops->evaluate)(field, points, datatype, B, D, H);
236: } else SETERRQ(PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "Not implemented for this type");
237: return 0;
238: }
240: /*@
241: DMFieldEvaluateFE - Evaluate the field and its derivatives on a set of points mapped from
242: quadrature points on a reference point. The derivatives are taken with respect to the
243: reference coordinates.
245: Not collective
247: Input Parameters:
248: + field - The `DMField` object
249: . cellIS - Index set for cells on which to evaluate the field
250: . points - The quadature containing the points in the reference cell at which to evaluate the field.
251: - datatype - The PetscDataType of the output arrays: either PETSC_REAL or PETSC_SCALAR.
252: If the field is complex and datatype is PETSC_REAL, the real part of the
253: field is returned.
255: Output Parameters:
256: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
257: If B is not NULL, the values of the field are written in this array, varying first by component,
258: then by point.
259: . D - pointer to data of size d * c * n * sizeof(datatype).
260: If D is not NULL, the values of the field's spatial derivatives are written in this array,
261: varying first by the partial derivative component, then by field component, then by point.
262: - H - pointer to data of size d * d * c * n * sizeof(datatype).
263: If H is not NULL, the values of the field's second spatial derivatives are written in this array,
264: varying first by the second partial derivative component, then by field component, then by point.
266: Level: intermediate
268: .seealso: `DMField`, `DM`, `DMFieldGetNumComponents()`, `DMFieldEvaluate()`, `DMFieldEvaluateFV()`
269: @*/
270: PetscErrorCode DMFieldEvaluateFE(DMField field, IS cellIS, PetscQuadrature points, PetscDataType datatype, void *B, void *D, void *H)
271: {
278: if (field->ops->evaluateFE) {
279: (*field->ops->evaluateFE)(field, cellIS, points, datatype, B, D, H);
280: } else SETERRQ(PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "Not implemented for this type");
281: return 0;
282: }
284: /*@
285: DMFieldEvaluateFV - Evaluate the mean of a field and its finite volume derivatives on a set of points.
287: Not collective
289: Input Parameters:
290: + field - The `DMField` object
291: . cellIS - Index set for cells on which to evaluate the field
292: - datatype - The PetscDataType of the output arrays: either `PETSC_REAL` or `PETSC_SCALAR`.
293: If the field is complex and datatype is `PETSC_REAL`, the real part of the
294: field is returned.
296: Output Parameters:
297: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
298: If B is not NULL, the values of the field are written in this array, varying first by component,
299: then by point.
300: . D - pointer to data of size d * c * n * sizeof(datatype).
301: If D is not NULL, the values of the field's spatial derivatives are written in this array,
302: varying first by the partial derivative component, then by field component, then by point.
303: - H - pointer to data of size d * d * c * n * sizeof(datatype).
304: If H is not NULL, the values of the field's second spatial derivatives are written in this array,
305: varying first by the second partial derivative component, then by field component, then by point.
307: Level: intermediate
309: .seealso: `DMField`, `IS`, `DMFieldGetNumComponents()`, `DMFieldEvaluate()`, `DMFieldEvaluateFE()`, `PetscDataType`
310: @*/
311: PetscErrorCode DMFieldEvaluateFV(DMField field, IS cellIS, PetscDataType datatype, void *B, void *D, void *H)
312: {
318: if (field->ops->evaluateFV) {
319: (*field->ops->evaluateFV)(field, cellIS, datatype, B, D, H);
320: } else SETERRQ(PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "Not implemented for this type");
321: return 0;
322: }
324: /*@
325: DMFieldGetDegree - Get the polynomial degree of a field when pulled back onto the
326: reference element
328: Not collective
330: Input Parameters:
331: + field - the `DMField` object
332: - cellIS - the index set of points over which we want know the invariance
334: Output Parameters:
335: + minDegree - the degree of the largest polynomial space contained in the field on each element
336: - maxDegree - the largest degree of the smallest polynomial space containing the field on any element
338: Level: intermediate
340: .seealso: `DMField`, `IS`, `DMFieldEvaluateFE()`
341: @*/
342: PetscErrorCode DMFieldGetDegree(DMField field, IS cellIS, PetscInt *minDegree, PetscInt *maxDegree)
343: {
349: if (minDegree) *minDegree = -1;
350: if (maxDegree) *maxDegree = PETSC_MAX_INT;
352: if (field->ops->getDegree) (*field->ops->getDegree)(field, cellIS, minDegree, maxDegree);
353: return 0;
354: }
356: /*@
357: DMFieldCreateDefaultQuadrature - Creates a quadrature sufficient to integrate the field on the selected
358: points via pullback onto the reference element
360: Not collective
362: Input Parameters:
363: + field - the `DMField` object
364: - pointIS - the index set of points over which we wish to integrate the field
366: Output Parameter:
367: . quad - a `PetscQuadrature` object
369: Level: developer
371: .seealso: `DMField`, `PetscQuadrature`, `IS`, `DMFieldEvaluteFE()`, `DMFieldGetDegree()`
372: @*/
373: PetscErrorCode DMFieldCreateDefaultQuadrature(DMField field, IS pointIS, PetscQuadrature *quad)
374: {
379: *quad = NULL;
380: PetscTryTypeMethod(field, createDefaultQuadrature, pointIS, quad);
381: return 0;
382: }
384: /*@C
385: DMFieldCreateFEGeom - Compute and create the geometric factors of a coordinate field
387: Not collective
389: Input Parameters:
390: + field - the `DMField` object
391: . pointIS - the index set of points over which we wish to integrate the field
392: . quad - the quadrature points at which to evaluate the geometric factors
393: - faceData - whether additional data for facets (the normal vectors and adjacent cells) should
394: be calculated
396: Output Parameter:
397: . geom - the geometric factors
399: Level: developer
401: .seealso: `DMField`, `PetscQuadrature`, `IS`, `PetscFEGeom`, `DMFieldEvaluateFE()`, `DMFieldCreateDefaulteQuadrature()`, `DMFieldGetDegree()`
402: @*/
403: PetscErrorCode DMFieldCreateFEGeom(DMField field, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
404: {
405: PetscInt dim, dE;
406: PetscInt nPoints;
407: PetscInt maxDegree;
408: PetscFEGeom *g;
413: ISGetLocalSize(pointIS, &nPoints);
414: dE = field->numComponents;
415: PetscFEGeomCreate(quad, nPoints, dE, faceData, &g);
416: DMFieldEvaluateFE(field, pointIS, quad, PETSC_REAL, g->v, g->J, NULL);
417: dim = g->dim;
418: if (dE > dim) {
419: /* space out J and make square Jacobians */
420: PetscInt i, j, k, N = g->numPoints * g->numCells;
422: for (i = N - 1; i >= 0; i--) {
423: PetscReal J[9];
425: for (j = 0; j < dE; j++) {
426: for (k = 0; k < dim; k++) J[j * dE + k] = g->J[i * dE * dim + j * dim + k];
427: }
428: switch (dim) {
429: case 0:
430: for (j = 0; j < dE; j++) {
431: for (k = 0; k < dE; k++) J[j * dE + k] = (j == k) ? 1. : 0.;
432: }
433: break;
434: case 1:
435: if (dE == 2) {
436: PetscReal norm = PetscSqrtReal(J[0] * J[0] + J[2] * J[2]);
438: J[1] = -J[2] / norm;
439: J[3] = J[0] / norm;
440: } else {
441: PetscReal inorm = 1. / PetscSqrtReal(J[0] * J[0] + J[3] * J[3] + J[6] * J[6]);
442: PetscReal x = J[0] * inorm;
443: PetscReal y = J[3] * inorm;
444: PetscReal z = J[6] * inorm;
446: if (x > 0.) {
447: PetscReal inv1pX = 1. / (1. + x);
449: J[1] = -y;
450: J[2] = -z;
451: J[4] = 1. - y * y * inv1pX;
452: J[5] = -y * z * inv1pX;
453: J[7] = -y * z * inv1pX;
454: J[8] = 1. - z * z * inv1pX;
455: } else {
456: PetscReal inv1mX = 1. / (1. - x);
458: J[1] = z;
459: J[2] = y;
460: J[4] = -y * z * inv1mX;
461: J[5] = 1. - y * y * inv1mX;
462: J[7] = 1. - z * z * inv1mX;
463: J[8] = -y * z * inv1mX;
464: }
465: }
466: break;
467: case 2: {
468: PetscReal inorm;
470: J[2] = J[3] * J[7] - J[6] * J[4];
471: J[5] = J[6] * J[1] - J[0] * J[7];
472: J[8] = J[0] * J[4] - J[3] * J[1];
474: inorm = 1. / PetscSqrtReal(J[2] * J[2] + J[5] * J[5] + J[8] * J[8]);
476: J[2] *= inorm;
477: J[5] *= inorm;
478: J[8] *= inorm;
479: } break;
480: }
481: for (j = 0; j < dE * dE; j++) g->J[i * dE * dE + j] = J[j];
482: }
483: }
484: PetscFEGeomComplete(g);
485: DMFieldGetDegree(field, pointIS, NULL, &maxDegree);
486: g->isAffine = (maxDegree <= 1) ? PETSC_TRUE : PETSC_FALSE;
487: if (faceData) (*field->ops->computeFaceData)(field, pointIS, quad, g);
488: *geom = g;
489: return 0;
490: }