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: /*@
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: /*@
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: PetscUseTypeMethod(field, evaluate, points, datatype, B, D, H);
243: PetscFunctionReturn(PETSC_SUCCESS);
244: }
246: /*@
247: DMFieldEvaluateFE - Evaluate the field and its derivatives on a set of points mapped from
248: quadrature points on a reference point. The derivatives are taken with respect to the
249: reference coordinates.
251: Not Collective
253: Input Parameters:
254: + field - The `DMField` object
255: . cellIS - Index set for cells on which to evaluate the field
256: . points - The quadature containing the points in the reference cell at which to evaluate the field.
257: - datatype - The PetscDataType of the output arrays: either `PETSC_REAL` or `PETSC_SCALAR`.
258: If the field is complex and datatype is `PETSC_REAL`, the real part of the
259: field is returned.
261: Output Parameters:
262: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
263: If B is not `NULL`, the values of the field are written in this array, varying first by component,
264: then by point.
265: . D - pointer to data of size d * c * n * sizeof(datatype).
266: If D is not `NULL`, the values of the field's spatial derivatives are written in this array,
267: varying first by the partial derivative component, then by field component, then by point.
268: - H - pointer to data of size d * d * c * n * sizeof(datatype).
269: If H is not `NULL`, the values of the field's second spatial derivatives are written in this array,
270: varying first by the second partial derivative component, then by field component, then by point.
272: Level: intermediate
274: .seealso: `DMField`, `DM`, `DMFieldGetNumComponents()`, `DMFieldEvaluate()`, `DMFieldEvaluateFV()`
275: @*/
276: PetscErrorCode DMFieldEvaluateFE(DMField field, IS cellIS, PetscQuadrature points, PetscDataType datatype, void *B, void *D, void *H)
277: {
278: PetscFunctionBegin;
282: if (B) PetscAssertPointer(B, 5);
283: if (D) PetscAssertPointer(D, 6);
284: if (H) PetscAssertPointer(H, 7);
285: PetscUseTypeMethod(field, evaluateFE, cellIS, points, datatype, B, D, H);
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
289: /*@
290: DMFieldEvaluateFV - Evaluate the mean of a field and its finite volume derivatives on a set of points.
292: Not Collective
294: Input Parameters:
295: + field - The `DMField` object
296: . cellIS - Index set for cells on which to evaluate the field
297: - datatype - The PetscDataType of the output arrays: either `PETSC_REAL` or `PETSC_SCALAR`.
298: If the field is complex and datatype is `PETSC_REAL`, the real part of the
299: field is returned.
301: Output Parameters:
302: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
303: If B is not `NULL`, the values of the field are written in this array, varying first by component,
304: then by point.
305: . D - pointer to data of size d * c * n * sizeof(datatype).
306: If D is not `NULL`, the values of the field's spatial derivatives are written in this array,
307: varying first by the partial derivative component, then by field component, then by point.
308: - H - pointer to data of size d * d * c * n * sizeof(datatype).
309: If H is not `NULL`, the values of the field's second spatial derivatives are written in this array,
310: varying first by the second partial derivative component, then by field component, then by point.
312: Level: intermediate
314: .seealso: `DMField`, `IS`, `DMFieldGetNumComponents()`, `DMFieldEvaluate()`, `DMFieldEvaluateFE()`, `PetscDataType`
315: @*/
316: PetscErrorCode DMFieldEvaluateFV(DMField field, IS cellIS, PetscDataType datatype, void *B, void *D, void *H)
317: {
318: PetscFunctionBegin;
321: if (B) PetscAssertPointer(B, 4);
322: if (D) PetscAssertPointer(D, 5);
323: if (H) PetscAssertPointer(H, 6);
324: PetscUseTypeMethod(field, evaluateFV, cellIS, datatype, B, D, H);
325: PetscFunctionReturn(PETSC_SUCCESS);
326: }
328: /*@
329: DMFieldGetDegree - Get the polynomial degree of a field when pulled back onto the
330: reference element
332: Not Collective
334: Input Parameters:
335: + field - the `DMField` object
336: - cellIS - the index set of points over which we want know the invariance
338: Output Parameters:
339: + minDegree - the degree of the largest polynomial space contained in the field on each element
340: - maxDegree - the largest degree of the smallest polynomial space containing the field on any element
342: Level: intermediate
344: .seealso: `DMField`, `IS`, `DMFieldEvaluateFE()`
345: @*/
346: PetscErrorCode DMFieldGetDegree(DMField field, IS cellIS, PetscInt *minDegree, PetscInt *maxDegree)
347: {
348: PetscFunctionBegin;
351: if (minDegree) PetscAssertPointer(minDegree, 3);
352: if (maxDegree) PetscAssertPointer(maxDegree, 4);
354: if (minDegree) *minDegree = -1;
355: if (maxDegree) *maxDegree = PETSC_MAX_INT;
357: PetscTryTypeMethod(field, getDegree, cellIS, minDegree, maxDegree);
358: PetscFunctionReturn(PETSC_SUCCESS);
359: }
361: /*@
362: DMFieldCreateDefaultQuadrature - Creates a quadrature sufficient to integrate the field on the selected
363: points via pullback onto the reference element
365: Not Collective
367: Input Parameters:
368: + field - the `DMField` object
369: - pointIS - the index set of points over which we wish to integrate the field
371: Output Parameter:
372: . quad - a `PetscQuadrature` object
374: Level: developer
376: .seealso: `DMField`, `PetscQuadrature`, `IS`, `DMFieldEvaluteFE()`, `DMFieldGetDegree()`
377: @*/
378: PetscErrorCode DMFieldCreateDefaultQuadrature(DMField field, IS pointIS, PetscQuadrature *quad)
379: {
380: PetscFunctionBegin;
383: PetscAssertPointer(quad, 3);
385: *quad = NULL;
386: PetscTryTypeMethod(field, createDefaultQuadrature, pointIS, quad);
387: PetscFunctionReturn(PETSC_SUCCESS);
388: }
390: /*@C
391: DMFieldCreateFEGeom - Compute and create the geometric factors of a coordinate field
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
398: . quad - the quadrature points at which to evaluate the geometric factors
399: - faceData - whether additional data for facets (the normal vectors and adjacent cells) should
400: be calculated
402: Output Parameter:
403: . geom - the geometric factors
405: Level: developer
407: .seealso: `DMField`, `PetscQuadrature`, `IS`, `PetscFEGeom`, `DMFieldEvaluateFE()`, `DMFieldCreateDefaulteQuadrature()`, `DMFieldGetDegree()`
408: @*/
409: PetscErrorCode DMFieldCreateFEGeom(DMField field, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
410: {
411: PetscInt dim, dE;
412: PetscInt nPoints;
413: PetscInt maxDegree;
414: PetscFEGeom *g;
416: PetscFunctionBegin;
420: PetscCall(ISGetLocalSize(pointIS, &nPoints));
421: dE = field->numComponents;
422: PetscCall(PetscFEGeomCreate(quad, nPoints, dE, faceData, &g));
423: PetscCall(DMFieldEvaluateFE(field, pointIS, quad, PETSC_REAL, g->v, g->J, NULL));
424: dim = g->dim;
425: if (dE > dim) {
426: /* space out J and make square Jacobians */
427: PetscInt i, j, k, N = g->numPoints * g->numCells;
429: for (i = N - 1; i >= 0; i--) {
430: PetscReal J[9] = {0};
432: for (j = 0; j < dE; j++) {
433: for (k = 0; k < dim; k++) J[j * dE + k] = g->J[i * dE * dim + j * dim + k];
434: }
435: switch (dim) {
436: case 0:
437: for (j = 0; j < dE; j++) {
438: for (k = 0; k < dE; k++) J[j * dE + k] = (j == k) ? 1. : 0.;
439: }
440: break;
441: case 1:
442: if (dE == 2) {
443: PetscReal norm = PetscSqrtReal(J[0] * J[0] + J[2] * J[2]);
445: J[1] = -J[2] / norm;
446: J[3] = J[0] / norm;
447: } else {
448: PetscReal inorm = 1. / PetscSqrtReal(J[0] * J[0] + J[3] * J[3] + J[6] * J[6]);
449: PetscReal x = J[0] * inorm;
450: PetscReal y = J[3] * inorm;
451: PetscReal z = J[6] * inorm;
453: if (x > 0.) {
454: PetscReal inv1pX = 1. / (1. + x);
456: J[1] = -y;
457: J[2] = -z;
458: J[4] = 1. - y * y * inv1pX;
459: J[5] = -y * z * inv1pX;
460: J[7] = -y * z * inv1pX;
461: J[8] = 1. - z * z * inv1pX;
462: } else {
463: PetscReal inv1mX = 1. / (1. - x);
465: J[1] = z;
466: J[2] = y;
467: J[4] = -y * z * inv1mX;
468: J[5] = 1. - y * y * inv1mX;
469: J[7] = 1. - z * z * inv1mX;
470: J[8] = -y * z * inv1mX;
471: }
472: }
473: break;
474: case 2: {
475: PetscReal inorm;
477: J[2] = J[3] * J[7] - J[6] * J[4];
478: J[5] = J[6] * J[1] - J[0] * J[7];
479: J[8] = J[0] * J[4] - J[3] * J[1];
481: inorm = 1. / PetscSqrtReal(J[2] * J[2] + J[5] * J[5] + J[8] * J[8]);
483: J[2] *= inorm;
484: J[5] *= inorm;
485: J[8] *= inorm;
486: } break;
487: }
488: for (j = 0; j < dE * dE; j++) g->J[i * dE * dE + j] = J[j];
489: }
490: }
491: PetscCall(PetscFEGeomComplete(g));
492: PetscCall(DMFieldGetDegree(field, pointIS, NULL, &maxDegree));
493: g->isAffine = (maxDegree <= 1) ? PETSC_TRUE : PETSC_FALSE;
494: if (faceData) PetscUseTypeMethod(field, computeFaceData, pointIS, quad, g);
495: *geom = g;
496: PetscFunctionReturn(PETSC_SUCCESS);
497: }