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: {
15: DMField b;
19: DMFieldInitializePackage();
21: PetscHeaderCreate(b,DMFIELD_CLASSID,"DMField","Field over DM","DM",PetscObjectComm((PetscObject)dm),DMFieldDestroy,DMFieldView);
22: PetscObjectReference((PetscObject)dm);
23: b->dm = dm;
24: b->continuity = continuity;
25: b->numComponents = numComponents;
26: *field = b;
27: return 0;
28: }
30: /*@
31: DMFieldDestroy - destroy a DMField
33: Collective
35: Input Parameter:
36: . field - address of DMField
38: Level: advanced
40: .seealso: DMFieldCreate()
41: @*/
42: PetscErrorCode DMFieldDestroy(DMField *field)
43: {
44: if (!*field) return 0;
46: if (--((PetscObject)(*field))->refct > 0) {*field = NULL; return 0;}
47: if ((*field)->ops->destroy) (*(*field)->ops->destroy)(*field);
48: DMDestroy(&((*field)->dm));
49: PetscHeaderDestroy(field);
50: return 0;
51: }
53: /*@C
54: DMFieldView - view a DMField
56: Collective
58: Input Parameters:
59: + field - DMField
60: - viewer - viewer to display field, for example PETSC_VIEWER_STDOUT_WORLD
62: Level: advanced
64: .seealso: DMFieldCreate()
65: @*/
66: PetscErrorCode DMFieldView(DMField field,PetscViewer viewer)
67: {
68: PetscBool iascii;
71: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)field),&viewer);
74: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
75: if (iascii) {
76: PetscObjectPrintClassNamePrefixType((PetscObject)field,viewer);
77: PetscViewerASCIIPushTab(viewer);
78: PetscViewerASCIIPrintf(viewer,"%D components\n",field->numComponents);
79: PetscViewerASCIIPrintf(viewer,"%s continuity\n",DMFieldContinuities[field->continuity]);
80: PetscViewerPushFormat(viewer,PETSC_VIEWER_DEFAULT);
81: DMView(field->dm,viewer);
82: PetscViewerPopFormat(viewer);
83: }
84: if (field->ops->view) (*field->ops->view)(field,viewer);
85: if (iascii) {
86: PetscViewerASCIIPopTab(viewer);
87: }
88: return 0;
89: }
91: /*@C
92: DMFieldSetType - set the DMField implementation
94: Collective on field
96: Input Parameters:
97: + field - the DMField context
98: - type - a known method
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: Level: advanced
108: .seealso: DMFieldType,
109: @*/
110: PetscErrorCode DMFieldSetType(DMField field,DMFieldType type)
111: {
112: PetscBool match;
113: PetscErrorCode (*r)(DMField);
118: PetscObjectTypeCompare((PetscObject)field,type,&match);
119: if (match) return 0;
121: PetscFunctionListFind(DMFieldList,type,&r);
123: /* Destroy the previous private DMField context */
124: if (field->ops->destroy) (*(field)->ops->destroy)(field);
126: PetscMemzero(field->ops,sizeof(*field->ops));
127: PetscObjectChangeTypeName((PetscObject)field,type);
128: field->ops->create = r;
129: (*r)(field);
130: return 0;
131: }
133: /*@C
134: DMFieldGetType - Gets the DMField type name (as a string) from the DMField.
136: Not Collective
138: Input Parameter:
139: . field - The DMField context
141: Output Parameter:
142: . type - The DMField type name
144: Level: advanced
146: .seealso: DMFieldSetType()
147: @*/
148: PetscErrorCode DMFieldGetType(DMField field, DMFieldType *type)
149: {
152: DMFieldRegisterAll();
153: *type = ((PetscObject)field)->type_name;
154: return 0;
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: DMFieldEvaluate()
171: @*/
172: PetscErrorCode DMFieldGetNumComponents(DMField field, PetscInt *nc)
173: {
176: *nc = field->numComponents;
177: return 0;
178: }
180: /*@
181: DMFieldGetDM - Returns the DM for the manifold over which the field is defined.
183: Not collective
185: Input Parameter:
186: . field - The DMField object
188: Output Parameter:
189: . dm - The DM object
191: Level: intermediate
193: .seealso: DMFieldEvaluate()
194: @*/
195: PetscErrorCode DMFieldGetDM(DMField field, DM *dm)
196: {
199: *dm = field->dm;
200: return 0;
201: }
203: /*@
204: DMFieldEvaluate - Evaluate the field and its derivatives on a set of points
206: Collective on points
208: Input Parameters:
209: + field - The DMField object
210: . points - The points at which to evaluate the field. Should have size d x n,
211: where d is the coordinate dimension of the manifold and n is the number
212: of points
213: - datatype - The PetscDataType of the output arrays: either PETSC_REAL or PETSC_SCALAR.
214: If the field is complex and datatype is PETSC_REAL, the real part of the
215: field is returned.
217: Output Parameters:
218: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
219: If B is not NULL, the values of the field are written in this array, varying first by component,
220: then by point.
221: . D - pointer to data of size d * c * n * sizeof(datatype).
222: If D is not NULL, the values of the field's spatial derivatives are written in this array,
223: varying first by the partial derivative component, then by field component, then by point.
224: - H - pointer to data of size d * d * c * n * sizeof(datatype).
225: If H is not NULL, the values of the field's second spatial derivatives are written in this array,
226: varying first by the second partial derivative component, then by field component, then by point.
228: Level: intermediate
230: .seealso: DMFieldGetDM(), DMFieldGetNumComponents(), DMFieldEvaluateFE(), DMFieldEvaluateFV()
231: @*/
232: PetscErrorCode DMFieldEvaluate(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
233: {
239: if (field->ops->evaluate) {
240: (*field->ops->evaluate) (field, points, datatype, B, D, H);
241: } else SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented for this type");
242: return 0;
243: }
245: /*@
246: DMFieldEvaluateFE - Evaluate the field and its derivatives on a set of points mapped from
247: quadrature points on a reference point. The derivatives are taken with respect to the
248: reference coordinates.
250: Not collective
252: Input Parameters:
253: + field - The DMField object
254: . cellIS - Index set for cells on which to evaluate the field
255: . points - The quadature containing the points in the reference cell at which to evaluate the field.
256: - datatype - The PetscDataType of the output arrays: either PETSC_REAL or PETSC_SCALAR.
257: If the field is complex and datatype is PETSC_REAL, the real part of the
258: field is returned.
260: Output Parameters:
261: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
262: If B is not NULL, the values of the field are written in this array, varying first by component,
263: then by point.
264: . D - pointer to data of size d * c * n * sizeof(datatype).
265: If D is not NULL, the values of the field's spatial derivatives are written in this array,
266: varying first by the partial derivative component, then by field component, then by point.
267: - H - pointer to data of size d * d * c * n * sizeof(datatype).
268: If H is not NULL, the values of the field's second spatial derivatives are written in this array,
269: varying first by the second partial derivative component, then by field component, then by point.
271: Level: intermediate
273: .seealso: DMFieldGetNumComponents(), DMFieldEvaluate(), DMFieldEvaluateFV()
274: @*/
275: PetscErrorCode DMFieldEvaluateFE(DMField field, IS cellIS, PetscQuadrature points, PetscDataType datatype, void *B, void *D, void *H)
276: {
283: if (field->ops->evaluateFE) {
284: (*field->ops->evaluateFE) (field, cellIS, points, datatype, B, D, H);
285: } else SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented for this type");
286: return 0;
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: DMFieldGetNumComponents(), DMFieldEvaluate(), DMFieldEvaluateFE()
315: @*/
316: PetscErrorCode DMFieldEvaluateFV(DMField field, IS cellIS, PetscDataType datatype, void *B, void *D, void *H)
317: {
323: if (field->ops->evaluateFV) {
324: (*field->ops->evaluateFV) (field, cellIS, datatype, B, D, H);
325: } else SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented for this type");
326: return 0;
327: }
329: /*@
330: DMFieldGetDegree - Get the polynomial degree of a field when pulled back onto the
331: reference element
333: Not collective
335: Input Parameters:
336: + field - the DMField object
337: - cellIS - the index set of points over which we want know the invariance
339: Output Parameters:
340: + minDegree - the degree of the largest polynomial space contained in the field on each element
341: - maxDegree - the largest degree of the smallest polynomial space containing the field on any element
343: Level: intermediate
345: .seealso: DMFieldEvaluateFE()
346: @*/
347: PetscErrorCode DMFieldGetDegree(DMField field, IS cellIS, PetscInt *minDegree, PetscInt *maxDegree)
348: {
354: if (minDegree) *minDegree = -1;
355: if (maxDegree) *maxDegree = PETSC_MAX_INT;
357: if (field->ops->getDegree) {
358: (*field->ops->getDegree) (field,cellIS,minDegree,maxDegree);
359: }
360: return 0;
361: }
363: /*@
364: DMFieldCreateDefaultQuadrature - Creates a quadrature sufficient to integrate the field on the selected
365: points via pullback onto the reference element
367: Not collective
369: Input Parameters:
370: + field - the DMField object
371: - pointIS - the index set of points over which we wish to integrate the field
373: Output Parameter:
374: . quad - a PetscQuadrature object
376: Level: developer
378: .seealso: DMFieldEvaluteFE(), DMFieldGetDegree()
379: @*/
380: PetscErrorCode DMFieldCreateDefaultQuadrature(DMField field, IS pointIS, PetscQuadrature *quad)
381: {
386: *quad = NULL;
387: if (field->ops->createDefaultQuadrature) {
388: (*field->ops->createDefaultQuadrature)(field, pointIS, quad);
389: }
390: return 0;
391: }
393: /*@C
394: DMFieldCreateFEGeom - Compute and create the geometric factors of a coordinate field
396: Not collective
398: Input Parameters:
399: + field - the DMField object
400: . pointIS - the index set of points over which we wish to integrate the field
401: . quad - the quadrature points at which to evaluate the geometric factors
402: - faceData - whether additional data for facets (the normal vectors and adjacent cells) should
403: be calculated
405: Output Parameter:
406: . geom - the geometric factors
408: Level: developer
410: .seealso: DMFieldEvaluateFE(), DMFieldCreateDefaulteQuadrature(), DMFieldGetDegree()
411: @*/
412: PetscErrorCode DMFieldCreateFEGeom(DMField field, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
413: {
414: PetscInt dim, dE;
415: PetscInt nPoints;
416: PetscInt maxDegree;
417: PetscFEGeom *g;
422: ISGetLocalSize(pointIS,&nPoints);
423: dE = field->numComponents;
424: PetscFEGeomCreate(quad,nPoints,dE,faceData,&g);
425: DMFieldEvaluateFE(field,pointIS,quad,PETSC_REAL,g->v,g->J,NULL);
426: dim = g->dim;
427: if (dE > dim) {
428: /* space out J and make square Jacobians */
429: PetscInt i, j, k, N = g->numPoints * g->numCells;
431: for (i = N-1; i >= 0; i--) {
432: PetscReal J[9];
434: for (j = 0; j < dE; j++) {
435: for (k = 0; k < dim; k++) {
436: J[j*dE + k] = g->J[i*dE*dim + j*dim + k];
437: }
438: }
439: switch (dim) {
440: case 0:
441: for (j = 0; j < dE; j++) {
442: for (k = 0; k < dE; k++) {
443: J[j * dE + k] = (j == k) ? 1. : 0.;
444: }
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; J[2] = -z;
463: J[4] = 1. - y*y*inv1pX; J[5] = -y*z*inv1pX;
464: J[7] = -y*z*inv1pX; J[8] = 1. - z*z*inv1pX;
465: } else {
466: PetscReal inv1mX = 1./ (1. - x);
468: J[1] = z; J[2] = y;
469: J[4] = -y*z*inv1mX; J[5] = 1. - y*y*inv1mX;
470: J[7] = 1. - z*z*inv1mX; J[8] = -y*z*inv1mX;
471: }
472: }
473: break;
474: case 2:
475: {
476: PetscReal inorm;
478: J[2] = J[3] * J[7] - J[6] * J[4];
479: J[5] = J[6] * J[1] - J[0] * J[7];
480: J[8] = J[0] * J[4] - J[3] * J[1];
482: inorm = 1./ PetscSqrtReal(J[2]*J[2] + J[5]*J[5] + J[8]*J[8]);
484: J[2] *= inorm;
485: J[5] *= inorm;
486: J[8] *= inorm;
487: }
488: break;
489: }
490: for (j = 0; j < dE*dE; j++) {
491: g->J[i*dE*dE + j] = J[j];
492: }
493: }
494: }
495: PetscFEGeomComplete(g);
496: DMFieldGetDegree(field,pointIS,NULL,&maxDegree);
497: g->isAffine = (maxDegree <= 1) ? PETSC_TRUE : PETSC_FALSE;
498: if (faceData) {
500: (*field->ops->computeFaceData) (field, pointIS, quad, g);
501: }
502: *geom = g;
503: return 0;
504: }