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: /*@
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 isascii;
69: PetscFunctionBegin;
71: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)field), &viewer));
73: PetscCheckSameComm(field, 1, viewer, 2);
74: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
75: if (isascii) {
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 (isascii) 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, see `DMFieldType`
98: Level: advanced
100: .seealso: `DMField`, `DMFieldGetType()`, `DMFieldType`
101: @*/
102: PetscErrorCode DMFieldSetType(DMField field, DMFieldType type)
103: {
104: PetscBool match;
105: PetscErrorCode (*r)(DMField);
107: PetscFunctionBegin;
109: PetscAssertPointer(type, 2);
111: PetscCall(PetscObjectTypeCompare((PetscObject)field, type, &match));
112: if (match) PetscFunctionReturn(PETSC_SUCCESS);
114: PetscCall(PetscFunctionListFind(DMFieldList, type, &r));
115: PetscCheck(r, PetscObjectComm((PetscObject)field), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unable to find requested DMField type %s", type);
116: /* Destroy the previous private DMField context */
117: PetscTryTypeMethod(field, destroy);
119: PetscCall(PetscMemzero(field->ops, sizeof(*field->ops)));
120: PetscCall(PetscObjectChangeTypeName((PetscObject)field, type));
121: field->ops->create = r;
122: PetscCall((*r)(field));
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: /*@
127: DMFieldGetType - Gets the `DMFieldType` name (as a string) from the `DMField`.
129: Not Collective
131: Input Parameter:
132: . field - The `DMField` context
134: Output Parameter:
135: . type - The `DMFieldType` name
137: Level: advanced
139: Note:
140: `type` should not be retained for later use as it will be an invalid pointer if the `DMFieldType` of `field` is changed.
142: .seealso: `DMField`, `DMFieldSetType()`, `DMFieldType`, `PetscObjectTypeCompare()`, `PetscObjectTypeCompareAny()`
143: @*/
144: PetscErrorCode DMFieldGetType(DMField field, DMFieldType *type)
145: {
146: PetscFunctionBegin;
148: PetscAssertPointer(type, 2);
149: PetscCall(DMFieldRegisterAll());
150: *type = ((PetscObject)field)->type_name;
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: /*@
155: DMFieldGetNumComponents - Returns the number of components in the field
157: Not Collective
159: Input Parameter:
160: . field - The `DMField` object
162: Output Parameter:
163: . nc - The number of field components
165: Level: intermediate
167: .seealso: `DMField`, `DMFieldEvaluate()`
168: @*/
169: PetscErrorCode DMFieldGetNumComponents(DMField field, PetscInt *nc)
170: {
171: PetscFunctionBegin;
173: PetscAssertPointer(nc, 2);
174: *nc = field->numComponents;
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: /*@
179: DMFieldGetDM - Returns the `DM` for the manifold over which the field is defined.
181: Not Collective
183: Input Parameter:
184: . field - The `DMField` object
186: Output Parameter:
187: . dm - The `DM` object
189: Level: intermediate
191: .seealso: `DMField`, `DM`, `DMFieldEvaluate()`
192: @*/
193: PetscErrorCode DMFieldGetDM(DMField field, DM *dm)
194: {
195: PetscFunctionBegin;
197: PetscAssertPointer(dm, 2);
198: *dm = field->dm;
199: PetscFunctionReturn(PETSC_SUCCESS);
200: }
202: /*@
203: DMFieldEvaluate - Evaluate the field and its derivatives on a set of points
205: Collective
207: Input Parameters:
208: + field - The `DMField` object
209: . points - The points at which to evaluate the field. Should have size d x n,
210: where d is the coordinate dimension of the manifold and n is the number
211: of points
212: - datatype - The PetscDataType of the output arrays: either `PETSC_REAL` or `PETSC_SCALAR`.
213: If the field is complex and datatype is `PETSC_REAL`, the real part of the
214: field is returned.
216: Output Parameters:
217: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
218: If B is not `NULL`, the values of the field are written in this array, varying first by component,
219: then by point.
220: . D - pointer to data of size d * c * n * sizeof(datatype).
221: If `D` is not `NULL`, the values of the field's spatial derivatives are written in this array,
222: varying first by the partial derivative component, then by field component, then by point.
223: - H - pointer to data of size d * d * c * n * sizeof(datatype).
224: If `H` is not `NULL`, the values of the field's second spatial derivatives are written in this array,
225: varying first by the second partial derivative component, then by field component, then by point.
227: Level: intermediate
229: .seealso: `DMField`, `DMFieldGetDM()`, `DMFieldGetNumComponents()`, `DMFieldEvaluateFE()`, `DMFieldEvaluateFV()`, `PetscDataType`
230: @*/
231: PetscErrorCode DMFieldEvaluate(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
232: {
233: PetscFunctionBegin;
236: if (B) PetscAssertPointer(B, 4);
237: if (D) PetscAssertPointer(D, 5);
238: if (H) PetscAssertPointer(H, 6);
239: PetscUseTypeMethod(field, evaluate, points, datatype, B, D, H);
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: /*@
244: DMFieldEvaluateFE - Evaluate the field and its derivatives on a set of points mapped from
245: quadrature points on a reference point. The derivatives are taken with respect to the
246: reference coordinates.
248: Not Collective
250: Input Parameters:
251: + field - The `DMField` object
252: . cellIS - Index set for cells on which to evaluate the field
253: . points - The quadature containing the points in the reference cell at which to evaluate the field.
254: - datatype - The PetscDataType of the output arrays: either `PETSC_REAL` or `PETSC_SCALAR`.
255: If the field is complex and datatype is `PETSC_REAL`, the real part of the
256: field is returned.
258: Output Parameters:
259: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
260: If B is not `NULL`, the values of the field are written in this array, varying first by component,
261: then by point.
262: . D - pointer to data of size d * c * n * sizeof(datatype).
263: If D is not `NULL`, the values of the field's spatial derivatives are written in this array,
264: varying first by the partial derivative component, then by field component, then by point.
265: - H - pointer to data of size d * d * c * n * sizeof(datatype).
266: If H is not `NULL`, the values of the field's second spatial derivatives are written in this array,
267: varying first by the second partial derivative component, then by field component, then by point.
269: Level: intermediate
271: .seealso: `DMField`, `DM`, `DMFieldGetNumComponents()`, `DMFieldEvaluate()`, `DMFieldEvaluateFV()`
272: @*/
273: PetscErrorCode DMFieldEvaluateFE(DMField field, IS cellIS, PetscQuadrature points, PetscDataType datatype, void *B, void *D, void *H)
274: {
275: PetscFunctionBegin;
279: if (B) PetscAssertPointer(B, 5);
280: if (D) PetscAssertPointer(D, 6);
281: if (H) PetscAssertPointer(H, 7);
282: PetscUseTypeMethod(field, evaluateFE, cellIS, points, datatype, B, D, H);
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: /*@
287: DMFieldEvaluateFV - Evaluate the mean of a field and its finite volume derivatives on a set of points.
289: Not Collective
291: Input Parameters:
292: + field - The `DMField` object
293: . cellIS - Index set for cells on which to evaluate the field
294: - datatype - The PetscDataType of the output arrays: either `PETSC_REAL` or `PETSC_SCALAR`.
295: If the field is complex and datatype is `PETSC_REAL`, the real part of the
296: field is returned.
298: Output Parameters:
299: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
300: If B is not `NULL`, the values of the field are written in this array, varying first by component,
301: then by point.
302: . D - pointer to data of size d * c * n * sizeof(datatype).
303: If D is not `NULL`, the values of the field's spatial derivatives are written in this array,
304: varying first by the partial derivative component, then by field component, then by point.
305: - H - pointer to data of size d * d * c * n * sizeof(datatype).
306: If H is not `NULL`, the values of the field's second spatial derivatives are written in this array,
307: varying first by the second partial derivative component, then by field component, then by point.
309: Level: intermediate
311: .seealso: `DMField`, `IS`, `DMFieldGetNumComponents()`, `DMFieldEvaluate()`, `DMFieldEvaluateFE()`, `PetscDataType`
312: @*/
313: PetscErrorCode DMFieldEvaluateFV(DMField field, IS cellIS, PetscDataType datatype, void *B, void *D, void *H)
314: {
315: PetscFunctionBegin;
318: if (B) PetscAssertPointer(B, 4);
319: if (D) PetscAssertPointer(D, 5);
320: if (H) PetscAssertPointer(H, 6);
321: PetscUseTypeMethod(field, evaluateFV, cellIS, datatype, B, D, H);
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: /*@
326: DMFieldGetDegree - Get the polynomial degree of a field when pulled back onto the
327: reference element
329: Not Collective
331: Input Parameters:
332: + field - the `DMField` object
333: - cellIS - the index set of points over which we want know the invariance
335: Output Parameters:
336: + minDegree - the degree of the largest polynomial space contained in the field on each element
337: - maxDegree - the largest degree of the smallest polynomial space containing the field on any element
339: Level: intermediate
341: .seealso: `DMField`, `IS`, `DMFieldEvaluateFE()`
342: @*/
343: PetscErrorCode DMFieldGetDegree(DMField field, IS cellIS, PeOp PetscInt *minDegree, PeOp PetscInt *maxDegree)
344: {
345: PetscFunctionBegin;
348: if (minDegree) PetscAssertPointer(minDegree, 3);
349: if (maxDegree) PetscAssertPointer(maxDegree, 4);
351: if (minDegree) *minDegree = -1;
352: if (maxDegree) *maxDegree = PETSC_INT_MAX;
354: PetscTryTypeMethod(field, getDegree, cellIS, minDegree, maxDegree);
355: PetscFunctionReturn(PETSC_SUCCESS);
356: }
358: /*@
359: DMFieldCreateDefaultQuadrature - Creates a quadrature sufficient to integrate the field on the selected
360: points via pullback onto the reference element
362: Not Collective
364: Input Parameters:
365: + field - the `DMField` object
366: - pointIS - the index set of points over which we wish to integrate the field
368: Output Parameter:
369: . quad - a `PetscQuadrature` object
371: Level: developer
373: .seealso: `DMFieldCreateDefaultFaceQuadrature()`, `DMField`, `PetscQuadrature`, `IS`, `DMFieldEvaluteFE()`, `DMFieldGetDegree()`
374: @*/
375: PetscErrorCode DMFieldCreateDefaultQuadrature(DMField field, IS pointIS, PetscQuadrature *quad)
376: {
377: PetscFunctionBegin;
380: PetscAssertPointer(quad, 3);
382: *quad = NULL;
383: PetscTryTypeMethod(field, createDefaultQuadrature, pointIS, quad);
384: PetscFunctionReturn(PETSC_SUCCESS);
385: }
387: /*@
388: DMFieldCreateDefaultFaceQuadrature - Creates a quadrature sufficient to integrate the field on all faces of the selected cells via pullback onto the reference element
390: Not Collective
392: Input Parameters:
393: + field - the `DMField` object
394: - pointIS - the index set of points over which we wish to integrate the field over faces
396: Output Parameter:
397: . quad - a `PetscQuadrature` object
399: Level: developer
401: .seealso: `DMFieldCreateDefaultQuadrature()`, `DMField`, `PetscQuadrature`, `IS`, `DMFieldEvaluteFE()`, `DMFieldGetDegree()`
402: @*/
403: PetscErrorCode DMFieldCreateDefaultFaceQuadrature(DMField field, IS pointIS, PetscQuadrature *quad)
404: {
405: PetscFunctionBegin;
408: PetscAssertPointer(quad, 3);
410: *quad = NULL;
411: PetscTryTypeMethod(field, createDefaultFaceQuadrature, pointIS, quad);
412: PetscFunctionReturn(PETSC_SUCCESS);
413: }
415: /*@C
416: DMFieldCreateFEGeom - Compute and create the geometric factors of a coordinate field
418: Not Collective
420: Input Parameters:
421: + field - the `DMField` object
422: . pointIS - the index set of points over which we wish to integrate the field
423: . quad - the quadrature points at which to evaluate the geometric factors
424: - mode - Type of geometry data to store
426: Output Parameter:
427: . geom - the geometric factors
429: Level: developer
431: Note:
432: For some modes, the normal vectors and adjacent cells are calculated
434: .seealso: `DMField`, `PetscQuadrature`, `IS`, `PetscFEGeom`, `DMFieldEvaluateFE()`, `DMFieldCreateDefaulteQuadrature()`, `DMFieldGetDegree()`
435: @*/
436: PetscErrorCode DMFieldCreateFEGeom(DMField field, IS pointIS, PetscQuadrature quad, PetscFEGeomMode mode, PetscFEGeom **geom)
437: {
438: PetscBool faceData = mode == PETSC_FEGEOM_BOUNDARY || mode == PETSC_FEGEOM_COHESIVE ? PETSC_TRUE : PETSC_FALSE;
439: PetscInt dim, dE;
440: PetscInt nPoints;
441: PetscInt maxDegree;
442: PetscFEGeom *g;
444: PetscFunctionBegin;
448: PetscCall(ISGetLocalSize(pointIS, &nPoints));
449: dE = field->numComponents;
450: PetscCall(PetscFEGeomCreate(quad, nPoints, dE, mode, &g));
451: PetscCall(DMFieldEvaluateFE(field, pointIS, quad, PETSC_REAL, g->v, g->J, NULL));
452: dim = g->dim;
453: if (dE > dim) {
454: /* space out J and make square Jacobians */
455: PetscInt i, j, k, N = g->numPoints * g->numCells;
457: for (i = N - 1; i >= 0; i--) {
458: PetscReal J[16] = {0};
460: for (j = 0; j < dE; j++) {
461: for (k = 0; k < dim; k++) J[j * dE + k] = g->J[i * dE * dim + j * dim + k];
462: }
463: switch (dim) {
464: case 0:
465: for (j = 0; j < dE; j++) {
466: for (k = 0; k < dE; k++) J[j * dE + k] = (j == k) ? 1. : 0.;
467: }
468: break;
469: case 1:
470: if (dE == 2) {
471: PetscReal norm = PetscSqrtReal(J[0] * J[0] + J[2] * J[2]);
473: J[1] = -J[2] / norm;
474: J[3] = J[0] / norm;
475: } else {
476: PetscReal inorm = 1. / PetscSqrtReal(J[0] * J[0] + J[3] * J[3] + J[6] * J[6]);
477: PetscReal x = J[0] * inorm;
478: PetscReal y = J[3] * inorm;
479: PetscReal z = J[6] * inorm;
481: if (x > 0.) {
482: PetscReal inv1pX = 1. / (1. + x);
484: J[1] = -y;
485: J[2] = -z;
486: J[4] = 1. - y * y * inv1pX;
487: J[5] = -y * z * inv1pX;
488: J[7] = -y * z * inv1pX;
489: J[8] = 1. - z * z * inv1pX;
490: } else {
491: PetscReal inv1mX = 1. / (1. - x);
493: J[1] = z;
494: J[2] = y;
495: J[4] = -y * z * inv1mX;
496: J[5] = 1. - y * y * inv1mX;
497: J[7] = 1. - z * z * inv1mX;
498: J[8] = -y * z * inv1mX;
499: }
500: }
501: break;
502: case 2: {
503: PetscReal inorm;
505: J[2] = J[3] * J[7] - J[6] * J[4];
506: J[5] = J[6] * J[1] - J[0] * J[7];
507: J[8] = J[0] * J[4] - J[3] * J[1];
509: inorm = 1. / PetscSqrtReal(J[2] * J[2] + J[5] * J[5] + J[8] * J[8]);
511: J[2] *= inorm;
512: J[5] *= inorm;
513: J[8] *= inorm;
514: } break;
515: }
516: for (j = 0; j < dE * dE; j++) g->J[i * dE * dE + j] = J[j];
517: }
518: }
519: PetscCall(PetscFEGeomComplete(g));
520: PetscCall(DMFieldGetDegree(field, pointIS, NULL, &maxDegree));
521: g->isAffine = (maxDegree <= 1) ? PETSC_TRUE : PETSC_FALSE;
522: if (faceData) PetscUseTypeMethod(field, computeFaceData, pointIS, quad, g);
523: *geom = g;
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }