Actual source code: dmfield.c
petsc-3.12.5 2020-03-29
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: 0
11: };
13: PETSC_INTERN PetscErrorCode DMFieldCreate(DM dm,PetscInt numComponents,DMFieldContinuity continuity,DMField *field)
14: {
16: DMField b;
21: DMFieldInitializePackage();
23: PetscHeaderCreate(b,DMFIELD_CLASSID,"DMField","Field over DM","DM",PetscObjectComm((PetscObject)dm),DMFieldDestroy,DMFieldView);
24: PetscObjectReference((PetscObject)dm);
25: b->dm = dm;
26: b->continuity = continuity;
27: b->numComponents = numComponents;
28: *field = b;
29: return(0);
30: }
32: /*@
33: DMFieldDestroy - destroy a DMField
35: Collective
37: Input Arguments:
38: . field - address of DMField
40: Level: advanced
42: .seealso: DMFieldCreate()
43: @*/
44: PetscErrorCode DMFieldDestroy(DMField *field)
45: {
49: if (!*field) return(0);
51: if (--((PetscObject)(*field))->refct > 0) {*field = 0; return(0);}
52: if ((*field)->ops->destroy) {(*(*field)->ops->destroy)(*field);}
53: DMDestroy(&((*field)->dm));
54: PetscHeaderDestroy(field);
55: return(0);
56: }
58: /*@C
59: DMFieldView - view a DMField
61: Collective
63: Input Arguments:
64: + field - DMField
65: - viewer - viewer to display field, for example PETSC_VIEWER_STDOUT_WORLD
67: Level: advanced
69: .seealso: DMFieldCreate()
70: @*/
71: PetscErrorCode DMFieldView(DMField field,PetscViewer viewer)
72: {
73: PetscErrorCode ierr;
74: PetscBool iascii;
78: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)field),&viewer);}
81: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
82: if (iascii) {
83: PetscObjectPrintClassNamePrefixType((PetscObject)field,viewer);
84: PetscViewerASCIIPushTab(viewer);
85: PetscViewerASCIIPrintf(viewer,"%D components\n",field->numComponents);
86: PetscViewerASCIIPrintf(viewer,"%s continuity\n",DMFieldContinuities[field->continuity]);
87: PetscViewerPushFormat(viewer,PETSC_VIEWER_DEFAULT);
88: DMView(field->dm,viewer);
89: PetscViewerPopFormat(viewer);
90: }
91: if (field->ops->view) {(*field->ops->view)(field,viewer);}
92: if (iascii) {
93: PetscViewerASCIIPopTab(viewer);
94: }
95: return(0);
96: }
98: /*@C
99: DMFieldSetType - set the DMField implementation
101: Collective on field
103: Input Parameters:
104: + field - the DMField context
105: - type - a known method
107: Notes:
108: See "include/petscvec.h" for available methods (for instance)
109: + DMFIELDDA - a field defined only by its values at the corners of a DMDA
110: . DMFIELDDS - a field defined by a discretization over a mesh set with DMSetField()
111: - DMFIELDSHELL - a field defined by arbitrary callbacks
113: Level: advanced
115: .seealso: DMFieldType,
116: @*/
117: PetscErrorCode DMFieldSetType(DMField field,DMFieldType type)
118: {
119: PetscErrorCode ierr,(*r)(DMField);
120: PetscBool match;
126: PetscObjectTypeCompare((PetscObject)field,type,&match);
127: if (match) return(0);
129: PetscFunctionListFind(DMFieldList,type,&r);
130: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested DMField type %s",type);
131: /* Destroy the previous private DMField context */
132: if (field->ops->destroy) {
133: (*(field)->ops->destroy)(field);
134: }
135: PetscMemzero(field->ops,sizeof(*field->ops));
136: PetscObjectChangeTypeName((PetscObject)field,type);
137: field->ops->create = r;
138: (*r)(field);
139: return(0);
140: }
142: /*@C
143: DMFieldGetType - Gets the DMField type name (as a string) from the DMField.
145: Not Collective
147: Input Parameter:
148: . field - The DMField context
150: Output Parameter:
151: . type - The DMField type name
153: Level: advanced
155: .seealso: DMFieldSetType()
156: @*/
157: PetscErrorCode DMFieldGetType(DMField field, DMFieldType *type)
158: {
164: DMFieldRegisterAll();
165: *type = ((PetscObject)field)->type_name;
166: return(0);
167: }
169: /*@
170: DMFieldGetNumComponents - Returns the number of components in the field
172: Not collective
174: Input Parameter:
175: . field - The DMField object
177: Output Parameter:
178: . nc - The number of field components
180: Level: intermediate
182: .seealso: DMFieldEvaluate()
183: @*/
184: PetscErrorCode DMFieldGetNumComponents(DMField field, PetscInt *nc)
185: {
189: *nc = field->numComponents;
190: return(0);
191: }
193: /*@
194: DMFieldGetDM - Returns the DM for the manifold over which the field is defined.
196: Not collective
198: Input Parameter:
199: . field - The DMField object
201: Output Parameter:
202: . dm - The DM object
204: Level: intermediate
206: .seealso: DMFieldEvaluate()
207: @*/
208: PetscErrorCode DMFieldGetDM(DMField field, DM *dm)
209: {
213: *dm = field->dm;
214: return(0);
215: }
217: /*@
218: DMFieldEvaluate - Evaluate the field and its derivatives on a set of points
220: Collective on points
222: Input Parameter:
223: + field - The DMField object
224: . points - The points at which to evaluate the field. Should have size d x n,
225: where d is the coordinate dimension of the manifold and n is the number
226: of points
227: - datatype - The PetscDataType of the output arrays: either PETSC_REAL or PETSC_SCALAR.
228: If the field is complex and datatype is PETSC_REAL, the real part of the
229: field is returned.
232: Output Parameter:
233: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
234: If B is not NULL, the values of the field are written in this array, varying first by component,
235: then by point.
236: . D - pointer to data of size d * c * n * sizeof(datatype).
237: If D is not NULL, the values of the field's spatial derivatives are written in this array,
238: varying first by the partial derivative component, then by field component, then by point.
239: - H - pointer to data of size d * d * c * n * sizeof(datatype).
240: If H is not NULL, the values of the field's second spatial derivatives are written in this array,
241: varying first by the second partial derivative component, then by field component, then by point.
243: Level: intermediate
245: .seealso: DMFieldGetDM(), DMFieldGetNumComponents(), DMFieldEvaluateFE(), DMFieldEvaluateFV()
246: @*/
247: PetscErrorCode DMFieldEvaluate(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
248: {
257: if (field->ops->evaluate) {
258: (*field->ops->evaluate) (field, points, datatype, B, D, H);
259: } else SETERRQ (PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented for this type");
260: return(0);
261: }
263: /*@
264: DMFieldEvaluateFE - Evaluate the field and its derivatives on a set of points mapped from
265: quadrature points on a reference point. The derivatives are taken with respect to the
266: reference coordinates.
268: Not collective
270: Input Parameter:
271: + field - The DMField object
272: . cellIS - Index set for cells on which to evaluate the field
273: . points - The quadature containing the points in the reference cell at which to evaluate the field.
274: - datatype - The PetscDataType of the output arrays: either PETSC_REAL or PETSC_SCALAR.
275: If the field is complex and datatype is PETSC_REAL, the real part of the
276: field is returned.
279: Output Parameter:
280: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
281: If B is not NULL, the values of the field are written in this array, varying first by component,
282: then by point.
283: . D - pointer to data of size d * c * n * sizeof(datatype).
284: If D is not NULL, the values of the field's spatial derivatives are written in this array,
285: varying first by the partial derivative component, then by field component, then by point.
286: - H - pointer to data of size d * d * c * n * sizeof(datatype).
287: If H is not NULL, the values of the field's second spatial derivatives are written in this array,
288: varying first by the second partial derivative component, then by field component, then by point.
290: Level: intermediate
292: .seealso: DMFieldGetNumComponents(), DMFieldEvaluate(), DMFieldEvaluateFV()
293: @*/
294: PetscErrorCode DMFieldEvaluateFE(DMField field, IS cellIS, PetscQuadrature points, PetscDataType datatype, void *B, void *D, void *H)
295: {
305: if (field->ops->evaluateFE) {
306: (*field->ops->evaluateFE) (field, cellIS, points, datatype, B, D, H);
307: } else SETERRQ (PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented for this type");
308: return(0);
309: }
311: /*@
312: DMFieldEvaluateFV - Evaluate the mean of a field and its finite volume derivatives on a set of points.
314: Not collective
316: Input Parameter:
317: + field - The DMField object
318: . cellIS - Index set for cells on which to evaluate the field
319: - datatype - The PetscDataType of the output arrays: either PETSC_REAL or PETSC_SCALAR.
320: If the field is complex and datatype is PETSC_REAL, the real part of the
321: field is returned.
324: Output Parameter:
325: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
326: If B is not NULL, the values of the field are written in this array, varying first by component,
327: then by point.
328: . D - pointer to data of size d * c * n * sizeof(datatype).
329: If D is not NULL, the values of the field's spatial derivatives are written in this array,
330: varying first by the partial derivative component, then by field component, then by point.
331: - H - pointer to data of size d * d * c * n * sizeof(datatype).
332: If H is not NULL, the values of the field's second spatial derivatives are written in this array,
333: varying first by the second partial derivative component, then by field component, then by point.
335: Level: intermediate
337: .seealso: DMFieldGetNumComponents(), DMFieldEvaluate(), DMFieldEvaluateFE()
338: @*/
339: PetscErrorCode DMFieldEvaluateFV(DMField field, IS cellIS, PetscDataType datatype, void *B, void *D, void *H)
340: {
349: if (field->ops->evaluateFV) {
350: (*field->ops->evaluateFV) (field, cellIS, datatype, B, D, H);
351: } else SETERRQ (PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented for this type");
352: return(0);
353: }
355: /*@
356: DMFieldGetDegree - Get the polynomial degree of a field when pulled back onto the
357: reference element
359: Not collective
361: Input Arguments:
362: + field - the DMField object
363: - cellIS - the index set of points over which we want know the invariance
365: Output Arguments:
366: + minDegree - the degree of the largest polynomial space contained in the field on each element
367: - maxDegree - the largest degree of the smallest polynomial space containing the field on any element
369: Level: intermediate
371: .seealso: DMFieldEvaluateFE()
372: @*/
373: PetscErrorCode DMFieldGetDegree(DMField field, IS cellIS, PetscInt *minDegree, PetscInt *maxDegree)
374: {
383: if (minDegree) *minDegree = -1;
384: if (maxDegree) *maxDegree = PETSC_MAX_INT;
386: if (field->ops->getDegree) {
387: (*field->ops->getDegree) (field,cellIS,minDegree,maxDegree);
388: }
389: return(0);
390: }
392: /*@
393: DMFieldCreateDefaultQuadrature - Creates a quadrature sufficient to integrate the field on the selected
394: points via pullback onto the reference element
396: Not collective
398: Input Arguments:
399: + field - the DMField object
400: - pointIS - the index set of points over which we wish to integrate the field
402: Output Arguments:
403: . quad - a PetscQuadrature object
405: Level: developer
407: .seealso: DMFieldEvaluteFE(), DMFieldGetDegree()
408: @*/
409: PetscErrorCode DMFieldCreateDefaultQuadrature(DMField field, IS pointIS, PetscQuadrature *quad)
410: {
418: *quad = NULL;
419: if (field->ops->createDefaultQuadrature) {
420: (*field->ops->createDefaultQuadrature)(field, pointIS, quad);
421: }
422: return(0);
423: }
425: /*@C
426: DMFieldCreateFEGeom - Compute and create the geometric factors of a coordinate field
428: Not collective
430: Input Arguments:
431: + field - the DMField object
432: . pointIS - the index set of points over which we wish to integrate the field
433: . quad - the quadrature points at which to evaluate the geometric factors
434: - faceData - whether additional data for facets (the normal vectors and adjacent cells) should
435: be calculated
437: Output Arguments:
438: . geom - the geometric factors
440: Level: developer
442: .seealso: DMFieldEvaluateFE(), DMFieldCreateDefaulteQuadrature(), DMFieldGetDegree()
443: @*/
444: PetscErrorCode DMFieldCreateFEGeom(DMField field, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
445: {
446: PetscInt dim, dE;
447: PetscInt nPoints;
448: PetscInt maxDegree;
449: PetscFEGeom *g;
456: ISGetLocalSize(pointIS,&nPoints);
457: dE = field->numComponents;
458: PetscFEGeomCreate(quad,nPoints,dE,faceData,&g);
459: DMFieldEvaluateFE(field,pointIS,quad,PETSC_REAL,g->v,g->J,NULL);
460: dim = g->dim;
461: if (dE > dim) {
462: /* space out J and make square Jacobians */
463: PetscInt i, j, k, N = g->numPoints * g->numCells;
465: for (i = N-1; i >= 0; i--) {
466: PetscReal J[9];
468: for (j = 0; j < dE; j++) {
469: for (k = 0; k < dim; k++) {
470: J[j*dE + k] = g->J[i*dE*dim + j*dim + k];
471: }
472: }
473: switch (dim) {
474: case 0:
475: for (j = 0; j < dE; j++) {
476: for (k = 0; k < dE; k++) {
477: J[j * dE + k] = (j == k) ? 1. : 0.;
478: }
479: }
480: break;
481: case 1:
482: if (dE == 2) {
483: PetscReal norm = PetscSqrtReal(J[0] * J[0] + J[2] * J[2]);
485: J[1] = -J[2] / norm;
486: J[3] = J[0] / norm;
487: } else {
488: PetscReal inorm = 1./PetscSqrtReal(J[0] * J[0] + J[3] * J[3] + J[6] * J[6]);
489: PetscReal x = J[0] * inorm;
490: PetscReal y = J[3] * inorm;
491: PetscReal z = J[6] * inorm;
493: if (x > 0.) {
494: PetscReal inv1pX = 1./ (1. + x);
496: J[1] = -y; J[2] = -z;
497: J[4] = 1. - y*y*inv1pX; J[5] = -y*z*inv1pX;
498: J[7] = -y*z*inv1pX; J[8] = 1. - z*z*inv1pX;
499: } else {
500: PetscReal inv1mX = 1./ (1. - x);
502: J[1] = z; J[2] = y;
503: J[4] = -y*z*inv1mX; J[5] = 1. - y*y*inv1mX;
504: J[7] = 1. - z*z*inv1mX; J[8] = -y*z*inv1mX;
505: }
506: }
507: break;
508: case 2:
509: {
510: PetscReal inorm;
512: J[2] = J[3] * J[7] - J[6] * J[4];
513: J[5] = J[6] * J[1] - J[0] * J[7];
514: J[8] = J[0] * J[4] - J[3] * J[1];
516: inorm = 1./ PetscSqrtReal(J[2]*J[2] + J[5]*J[5] + J[8]*J[8]);
518: J[2] *= inorm;
519: J[5] *= inorm;
520: J[8] *= inorm;
521: }
522: break;
523: }
524: for (j = 0; j < dE*dE; j++) {
525: g->J[i*dE*dE + j] = J[j];
526: }
527: }
528: }
529: PetscFEGeomComplete(g);
530: DMFieldGetDegree(field,pointIS,NULL,&maxDegree);
531: g->isAffine = (maxDegree <= 1) ? PETSC_TRUE : PETSC_FALSE;
532: if (faceData) {
533: if (!field->ops->computeFaceData) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "DMField implementation does not compute face data\n");
534: (*field->ops->computeFaceData) (field, pointIS, quad, g);
535: }
536: *geom = g;
537: return(0);
538: }