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;
11: PetscFunctionBegin;
13: PetscAssertPointer(field, 4);
14: PetscCall(DMFieldInitializePackage());
16: PetscCall(PetscHeaderCreate(b, DMFIELD_CLASSID, "DMField", "Field over DM", "DM", PetscObjectComm((PetscObject)dm), DMFieldDestroy, DMFieldView));
17: PetscCall(PetscObjectReference((PetscObject)dm));
18: b->dm = dm;
19: b->continuity = continuity;
20: b->numComponents = numComponents;
21: *field = b;
22: PetscFunctionReturn(PETSC_SUCCESS);
23: }
25: /*@
26: DMFieldDestroy - destroy a `DMField`
28: Collective
30: Input Parameter:
31: . field - address of `DMField`
33: Level: advanced
35: .seealso: `DMField`, `DMFieldCreate()`
36: @*/
37: PetscErrorCode DMFieldDestroy(DMField *field)
38: {
39: PetscFunctionBegin;
40: if (!*field) PetscFunctionReturn(PETSC_SUCCESS);
42: if (--((PetscObject)(*field))->refct > 0) {
43: *field = NULL;
44: PetscFunctionReturn(PETSC_SUCCESS);
45: }
46: PetscTryTypeMethod((*field), destroy);
47: PetscCall(DMDestroy(&((*field)->dm)));
48: PetscCall(PetscHeaderDestroy(field));
49: PetscFunctionReturn(PETSC_SUCCESS);
50: }
52: /*@C
53: DMFieldView - view a `DMField`
55: Collective
57: Input Parameters:
58: + field - `DMField`
59: - viewer - viewer to display field, for example `PETSC_VIEWER_STDOUT_WORLD`
61: Level: advanced
63: .seealso: `DMField`, `DMFieldCreate()`
64: @*/
65: PetscErrorCode DMFieldView(DMField field, PetscViewer viewer)
66: {
67: PetscBool iascii;
69: PetscFunctionBegin;
71: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)field), &viewer));
73: PetscCheckSameComm(field, 1, viewer, 2);
74: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
75: if (iascii) {
76: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)field, viewer));
77: PetscCall(PetscViewerASCIIPushTab(viewer));
78: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " components\n", field->numComponents));
79: PetscCall(PetscViewerASCIIPrintf(viewer, "%s continuity\n", DMFieldContinuities[field->continuity]));
80: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_DEFAULT));
81: PetscCall(DMView(field->dm, viewer));
82: PetscCall(PetscViewerPopFormat(viewer));
83: }
84: PetscTryTypeMethod(field, view, viewer);
85: if (iascii) PetscCall(PetscViewerASCIIPopTab(viewer));
86: PetscFunctionReturn(PETSC_SUCCESS);
87: }
89: /*@C
90: DMFieldSetType - set the `DMField` implementation
92: Collective
94: Input Parameters:
95: + field - the `DMField` context
96: - type - a known method
98: Level: advanced
100: Notes:
101: See "include/petscvec.h" for available methods (for instance)
102: + `DMFIELDDA` - a field defined only by its values at the corners of a `DMDA`
103: . `DMFIELDDS` - a field defined by a discretization over a mesh set with `DMSetField()`
104: - `DMFIELDSHELL` - a field defined by arbitrary callbacks
106: .seealso: `DMField`, `DMFieldGetType()`, `DMFieldType`,
107: @*/
108: PetscErrorCode DMFieldSetType(DMField field, DMFieldType type)
109: {
110: PetscBool match;
111: PetscErrorCode (*r)(DMField);
113: PetscFunctionBegin;
115: PetscAssertPointer(type, 2);
117: PetscCall(PetscObjectTypeCompare((PetscObject)field, type, &match));
118: if (match) PetscFunctionReturn(PETSC_SUCCESS);
120: PetscCall(PetscFunctionListFind(DMFieldList, type, &r));
121: PetscCheck(r, PetscObjectComm((PetscObject)field), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unable to find requested DMField type %s", type);
122: /* Destroy the previous private DMField context */
123: PetscTryTypeMethod(field, destroy);
125: PetscCall(PetscMemzero(field->ops, sizeof(*field->ops)));
126: PetscCall(PetscObjectChangeTypeName((PetscObject)field, type));
127: field->ops->create = r;
128: PetscCall((*r)(field));
129: PetscFunctionReturn(PETSC_SUCCESS);
130: }
132: /*@C
133: DMFieldGetType - Gets the `DMFieldType` name (as a string) from the `DMField`.
135: Not Collective
137: Input Parameter:
138: . field - The `DMField` context
140: Output Parameter:
141: . type - The `DMFieldType` name
143: Level: advanced
145: .seealso: `DMField`, `DMFieldSetType()`, `DMFieldType`
146: @*/
147: PetscErrorCode DMFieldGetType(DMField field, DMFieldType *type)
148: {
149: PetscFunctionBegin;
151: PetscAssertPointer(type, 2);
152: PetscCall(DMFieldRegisterAll());
153: *type = ((PetscObject)field)->type_name;
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
157: /*@
158: DMFieldGetNumComponents - Returns the number of components in the field
160: Not Collective
162: Input Parameter:
163: . field - The `DMField` object
165: Output Parameter:
166: . nc - The number of field components
168: Level: intermediate
170: .seealso: `DMField`, `DMFieldEvaluate()`
171: @*/
172: PetscErrorCode DMFieldGetNumComponents(DMField field, PetscInt *nc)
173: {
174: PetscFunctionBegin;
176: PetscAssertPointer(nc, 2);
177: *nc = field->numComponents;
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }
181: /*@
182: DMFieldGetDM - Returns the `DM` for the manifold over which the field is defined.
184: Not Collective
186: Input Parameter:
187: . field - The `DMField` object
189: Output Parameter:
190: . dm - The `DM` object
192: Level: intermediate
194: .seealso: `DMField`, `DM`, `DMFieldEvaluate()`
195: @*/
196: PetscErrorCode DMFieldGetDM(DMField field, DM *dm)
197: {
198: PetscFunctionBegin;
200: PetscAssertPointer(dm, 2);
201: *dm = field->dm;
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: /*@
206: DMFieldEvaluate - Evaluate the field and its derivatives on a set of points
208: Collective
210: Input Parameters:
211: + field - The `DMField` object
212: . points - The points at which to evaluate the field. Should have size d x n,
213: where d is the coordinate dimension of the manifold and n is the number
214: of points
215: - datatype - The PetscDataType of the output arrays: either `PETSC_REAL` or `PETSC_SCALAR`.
216: If the field is complex and datatype is `PETSC_REAL`, the real part of the
217: field is returned.
219: Output Parameters:
220: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
221: If B is not NULL, the values of the field are written in this array, varying first by component,
222: then by point.
223: . D - pointer to data of size d * c * n * sizeof(datatype).
224: If D is not NULL, the values of the field's spatial derivatives are written in this array,
225: varying first by the partial derivative component, then by field component, then by point.
226: - H - pointer to data of size d * d * c * n * sizeof(datatype).
227: If H is not NULL, the values of the field's second spatial derivatives are written in this array,
228: varying first by the second partial derivative component, then by field component, then by point.
230: Level: intermediate
232: .seealso: `DMField`, `DMFieldGetDM()`, `DMFieldGetNumComponents()`, `DMFieldEvaluateFE()`, `DMFieldEvaluateFV()`, `PetscDataType`
233: @*/
234: PetscErrorCode DMFieldEvaluate(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
235: {
236: PetscFunctionBegin;
239: if (B) PetscAssertPointer(B, 4);
240: if (D) PetscAssertPointer(D, 5);
241: if (H) PetscAssertPointer(H, 6);
242: if (field->ops->evaluate) {
243: PetscCall((*field->ops->evaluate)(field, points, datatype, B, D, H));
244: } else SETERRQ(PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "Not implemented for this type");
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: /*@
249: DMFieldEvaluateFE - Evaluate the field and its derivatives on a set of points mapped from
250: quadrature points on a reference point. The derivatives are taken with respect to the
251: reference coordinates.
253: Not Collective
255: Input Parameters:
256: + field - The `DMField` object
257: . cellIS - Index set for cells on which to evaluate the field
258: . points - The quadature containing the points in the reference cell at which to evaluate the field.
259: - datatype - The PetscDataType of the output arrays: either `PETSC_REAL` or `PETSC_SCALAR`.
260: If the field is complex and datatype is `PETSC_REAL`, the real part of the
261: field is returned.
263: Output Parameters:
264: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
265: If B is not `NULL`, the values of the field are written in this array, varying first by component,
266: then by point.
267: . D - pointer to data of size d * c * n * sizeof(datatype).
268: If D is not `NULL`, the values of the field's spatial derivatives are written in this array,
269: varying first by the partial derivative component, then by field component, then by point.
270: - H - pointer to data of size d * d * c * n * sizeof(datatype).
271: If H is not `NULL`, the values of the field's second spatial derivatives are written in this array,
272: varying first by the second partial derivative component, then by field component, then by point.
274: Level: intermediate
276: .seealso: `DMField`, `DM`, `DMFieldGetNumComponents()`, `DMFieldEvaluate()`, `DMFieldEvaluateFV()`
277: @*/
278: PetscErrorCode DMFieldEvaluateFE(DMField field, IS cellIS, PetscQuadrature points, PetscDataType datatype, void *B, void *D, void *H)
279: {
280: PetscFunctionBegin;
284: if (B) PetscAssertPointer(B, 5);
285: if (D) PetscAssertPointer(D, 6);
286: if (H) PetscAssertPointer(H, 7);
287: if (field->ops->evaluateFE) {
288: PetscCall((*field->ops->evaluateFE)(field, cellIS, points, datatype, B, D, H));
289: } else SETERRQ(PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "Not implemented for this type");
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: /*@
294: DMFieldEvaluateFV - Evaluate the mean of a field and its finite volume derivatives on a set of points.
296: Not Collective
298: Input Parameters:
299: + field - The `DMField` object
300: . cellIS - Index set for cells on which to evaluate the field
301: - datatype - The PetscDataType of the output arrays: either `PETSC_REAL` or `PETSC_SCALAR`.
302: If the field is complex and datatype is `PETSC_REAL`, the real part of the
303: field is returned.
305: Output Parameters:
306: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
307: If B is not `NULL`, the values of the field are written in this array, varying first by component,
308: then by point.
309: . D - pointer to data of size d * c * n * sizeof(datatype).
310: If D is not `NULL`, the values of the field's spatial derivatives are written in this array,
311: varying first by the partial derivative component, then by field component, then by point.
312: - H - pointer to data of size d * d * c * n * sizeof(datatype).
313: If H is not `NULL`, the values of the field's second spatial derivatives are written in this array,
314: varying first by the second partial derivative component, then by field component, then by point.
316: Level: intermediate
318: .seealso: `DMField`, `IS`, `DMFieldGetNumComponents()`, `DMFieldEvaluate()`, `DMFieldEvaluateFE()`, `PetscDataType`
319: @*/
320: PetscErrorCode DMFieldEvaluateFV(DMField field, IS cellIS, PetscDataType datatype, void *B, void *D, void *H)
321: {
322: PetscFunctionBegin;
325: if (B) PetscAssertPointer(B, 4);
326: if (D) PetscAssertPointer(D, 5);
327: if (H) PetscAssertPointer(H, 6);
328: if (field->ops->evaluateFV) {
329: PetscCall((*field->ops->evaluateFV)(field, cellIS, datatype, B, D, H));
330: } else SETERRQ(PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "Not implemented for this type");
331: PetscFunctionReturn(PETSC_SUCCESS);
332: }
334: /*@
335: DMFieldGetDegree - Get the polynomial degree of a field when pulled back onto the
336: reference element
338: Not Collective
340: Input Parameters:
341: + field - the `DMField` object
342: - cellIS - the index set of points over which we want know the invariance
344: Output Parameters:
345: + minDegree - the degree of the largest polynomial space contained in the field on each element
346: - maxDegree - the largest degree of the smallest polynomial space containing the field on any element
348: Level: intermediate
350: .seealso: `DMField`, `IS`, `DMFieldEvaluateFE()`
351: @*/
352: PetscErrorCode DMFieldGetDegree(DMField field, IS cellIS, PetscInt *minDegree, PetscInt *maxDegree)
353: {
354: PetscFunctionBegin;
357: if (minDegree) PetscAssertPointer(minDegree, 3);
358: if (maxDegree) PetscAssertPointer(maxDegree, 4);
360: if (minDegree) *minDegree = -1;
361: if (maxDegree) *maxDegree = PETSC_MAX_INT;
363: if (field->ops->getDegree) PetscCall((*field->ops->getDegree)(field, cellIS, minDegree, maxDegree));
364: PetscFunctionReturn(PETSC_SUCCESS);
365: }
367: /*@
368: DMFieldCreateDefaultQuadrature - Creates a quadrature sufficient to integrate the field on the selected
369: points via pullback onto the reference element
371: Not Collective
373: Input Parameters:
374: + field - the `DMField` object
375: - pointIS - the index set of points over which we wish to integrate the field
377: Output Parameter:
378: . quad - a `PetscQuadrature` object
380: Level: developer
382: .seealso: `DMField`, `PetscQuadrature`, `IS`, `DMFieldEvaluteFE()`, `DMFieldGetDegree()`
383: @*/
384: PetscErrorCode DMFieldCreateDefaultQuadrature(DMField field, IS pointIS, PetscQuadrature *quad)
385: {
386: PetscFunctionBegin;
389: PetscAssertPointer(quad, 3);
391: *quad = NULL;
392: PetscTryTypeMethod(field, createDefaultQuadrature, pointIS, quad);
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }
396: /*@C
397: DMFieldCreateFEGeom - Compute and create the geometric factors of a coordinate field
399: Not Collective
401: Input Parameters:
402: + field - the `DMField` object
403: . pointIS - the index set of points over which we wish to integrate the field
404: . quad - the quadrature points at which to evaluate the geometric factors
405: - faceData - whether additional data for facets (the normal vectors and adjacent cells) should
406: be calculated
408: Output Parameter:
409: . geom - the geometric factors
411: Level: developer
413: .seealso: `DMField`, `PetscQuadrature`, `IS`, `PetscFEGeom`, `DMFieldEvaluateFE()`, `DMFieldCreateDefaulteQuadrature()`, `DMFieldGetDegree()`
414: @*/
415: PetscErrorCode DMFieldCreateFEGeom(DMField field, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
416: {
417: PetscInt dim, dE;
418: PetscInt nPoints;
419: PetscInt maxDegree;
420: PetscFEGeom *g;
422: PetscFunctionBegin;
426: PetscCall(ISGetLocalSize(pointIS, &nPoints));
427: dE = field->numComponents;
428: PetscCall(PetscFEGeomCreate(quad, nPoints, dE, faceData, &g));
429: PetscCall(DMFieldEvaluateFE(field, pointIS, quad, PETSC_REAL, g->v, g->J, NULL));
430: dim = g->dim;
431: if (dE > dim) {
432: /* space out J and make square Jacobians */
433: PetscInt i, j, k, N = g->numPoints * g->numCells;
435: for (i = N - 1; i >= 0; i--) {
436: PetscReal J[9] = {0};
438: for (j = 0; j < dE; j++) {
439: for (k = 0; k < dim; k++) J[j * dE + k] = g->J[i * dE * dim + j * dim + k];
440: }
441: switch (dim) {
442: case 0:
443: for (j = 0; j < dE; j++) {
444: for (k = 0; k < dE; k++) J[j * dE + k] = (j == k) ? 1. : 0.;
445: }
446: break;
447: case 1:
448: if (dE == 2) {
449: PetscReal norm = PetscSqrtReal(J[0] * J[0] + J[2] * J[2]);
451: J[1] = -J[2] / norm;
452: J[3] = J[0] / norm;
453: } else {
454: PetscReal inorm = 1. / PetscSqrtReal(J[0] * J[0] + J[3] * J[3] + J[6] * J[6]);
455: PetscReal x = J[0] * inorm;
456: PetscReal y = J[3] * inorm;
457: PetscReal z = J[6] * inorm;
459: if (x > 0.) {
460: PetscReal inv1pX = 1. / (1. + x);
462: J[1] = -y;
463: J[2] = -z;
464: J[4] = 1. - y * y * inv1pX;
465: J[5] = -y * z * inv1pX;
466: J[7] = -y * z * inv1pX;
467: J[8] = 1. - z * z * inv1pX;
468: } else {
469: PetscReal inv1mX = 1. / (1. - x);
471: J[1] = z;
472: J[2] = y;
473: J[4] = -y * z * inv1mX;
474: J[5] = 1. - y * y * inv1mX;
475: J[7] = 1. - z * z * inv1mX;
476: J[8] = -y * z * inv1mX;
477: }
478: }
479: break;
480: case 2: {
481: PetscReal inorm;
483: J[2] = J[3] * J[7] - J[6] * J[4];
484: J[5] = J[6] * J[1] - J[0] * J[7];
485: J[8] = J[0] * J[4] - J[3] * J[1];
487: inorm = 1. / PetscSqrtReal(J[2] * J[2] + J[5] * J[5] + J[8] * J[8]);
489: J[2] *= inorm;
490: J[5] *= inorm;
491: J[8] *= inorm;
492: } break;
493: }
494: for (j = 0; j < dE * dE; j++) g->J[i * dE * dE + j] = J[j];
495: }
496: }
497: PetscCall(PetscFEGeomComplete(g));
498: PetscCall(DMFieldGetDegree(field, pointIS, NULL, &maxDegree));
499: g->isAffine = (maxDegree <= 1) ? PETSC_TRUE : PETSC_FALSE;
500: if (faceData) PetscCall((*field->ops->computeFaceData)(field, pointIS, quad, g));
501: *geom = g;
502: PetscFunctionReturn(PETSC_SUCCESS);
503: }