Actual source code: dmfield.c
petsc-3.10.5 2019-03-28
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 DMField
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: .keywords: DMField, set, type
117: .seealso: DMFieldType,
118: @*/
119: PetscErrorCode DMFieldSetType(DMField field,DMFieldType type)
120: {
121: PetscErrorCode ierr,(*r)(DMField);
122: PetscBool match;
128: PetscObjectTypeCompare((PetscObject)field,type,&match);
129: if (match) return(0);
131: PetscFunctionListFind(DMFieldList,type,&r);
132: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested DMField type %s",type);
133: /* Destroy the previous private DMField context */
134: if (field->ops->destroy) {
135: (*(field)->ops->destroy)(field);
136: }
137: PetscMemzero(field->ops,sizeof(*field->ops));
138: PetscObjectChangeTypeName((PetscObject)field,type);
139: field->ops->create = r;
140: (*r)(field);
141: return(0);
142: }
144: /*@C
145: DMFieldGetType - Gets the DMField type name (as a string) from the DMField.
147: Not Collective
149: Input Parameter:
150: . field - The DMField context
152: Output Parameter:
153: . type - The DMField type name
155: Level: advanced
157: .keywords: DMField, get, type, name
158: .seealso: DMFieldSetType()
159: @*/
160: PetscErrorCode DMFieldGetType(DMField field, DMFieldType *type)
161: {
167: DMFieldRegisterAll();
168: *type = ((PetscObject)field)->type_name;
169: return(0);
170: }
172: /*@
173: DMFieldGetNumComponents - Returns the number of components in the field
175: Not collective
177: Input Parameter:
178: . field - The DMField object
180: Output Parameter:
181: . nc - The number of field components
183: Level: intermediate
185: .seealso: DMFieldEvaluate()
186: @*/
187: PetscErrorCode DMFieldGetNumComponents(DMField field, PetscInt *nc)
188: {
192: *nc = field->numComponents;
193: return(0);
194: }
196: /*@
197: DMFieldGetDM - Returns the DM for the manifold over which the field is defined.
199: Not collective
201: Input Parameter:
202: . field - The DMField object
204: Output Parameter:
205: . dm - The DM object
207: Level: intermediate
209: .seealso: DMFieldEvaluate()
210: @*/
211: PetscErrorCode DMFieldGetDM(DMField field, DM *dm)
212: {
216: *dm = field->dm;
217: return(0);
218: }
220: /*@
221: DMFieldEvaluate - Evaluate the field and its derivatives on a set of points
223: Collective on Vec
225: Input Parameter:
226: + field - The DMField object
227: . points - The points at which to evaluate the field. Should have size d x n,
228: where d is the coordinate dimension of the manifold and n is the number
229: of points
230: - datatype - The PetscDataType of the output arrays: either PETSC_REAL or PETSC_SCALAR.
231: If the field is complex and datatype is PETSC_REAL, the real part of the
232: field is returned.
235: Output Parameter:
236: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
237: If B is not NULL, the values of the field are written in this array, varying first by component,
238: then by point.
239: . D - pointer to data of size d * c * n * sizeof(datatype).
240: If D is not NULL, the values of the field's spatial derivatives are written in this array,
241: varying first by the partial derivative component, then by field component, then by point.
242: - H - pointer to data of size d * d * c * n * sizeof(datatype).
243: If H is not NULL, the values of the field's second spatial derivatives are written in this array,
244: varying first by the second partial derivative component, then by field component, then by point.
246: Level: intermediate
248: .seealso: DMFieldGetDM(), DMFieldGetNumComponents(), DMFieldEvaluateFE(), DMFieldEvaluateFV()
249: @*/
250: PetscErrorCode DMFieldEvaluate(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
251: {
260: if (field->ops->evaluate) {
261: (*field->ops->evaluate) (field, points, datatype, B, D, H);
262: } else SETERRQ (PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented for this type");
263: return(0);
264: }
266: /*@
267: DMFieldEvaluateFE - Evaluate the field and its derivatives on a set of points mapped from
268: quadrature points on a reference point. The derivatives are taken with respect to the
269: reference coordinates.
271: Not collective
273: Input Parameter:
274: + field - The DMField object
275: . cellIS - Index set for cells on which to evaluate the field
276: . points - The quadature containing the points in the reference cell at which to evaluate the field.
277: - datatype - The PetscDataType of the output arrays: either PETSC_REAL or PETSC_SCALAR.
278: If the field is complex and datatype is PETSC_REAL, the real part of the
279: field is returned.
282: Output Parameter:
283: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
284: If B is not NULL, the values of the field are written in this array, varying first by component,
285: then by point.
286: . D - pointer to data of size d * c * n * sizeof(datatype).
287: If D is not NULL, the values of the field's spatial derivatives are written in this array,
288: varying first by the partial derivative component, then by field component, then by point.
289: - H - pointer to data of size d * d * c * n * sizeof(datatype).
290: If H is not NULL, the values of the field's second spatial derivatives are written in this array,
291: varying first by the second partial derivative component, then by field component, then by point.
293: Level: intermediate
295: .seealso: DMFieldGetNumComponents(), DMFieldEvaluate(), DMFieldEvaluateFV()
296: @*/
297: PetscErrorCode DMFieldEvaluateFE(DMField field, IS cellIS, PetscQuadrature points, PetscDataType datatype, void *B, void *D, void *H)
298: {
308: if (field->ops->evaluateFE) {
309: (*field->ops->evaluateFE) (field, cellIS, points, datatype, B, D, H);
310: } else SETERRQ (PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented for this type");
311: return(0);
312: }
314: /*@
315: DMFieldEvaluateFV - Evaluate the mean of a field and its finite volume derivatives on a set of points.
317: Not collective
319: Input Parameter:
320: + field - The DMField object
321: . cellIS - Index set for cells on which to evaluate the field
322: - datatype - The PetscDataType of the output arrays: either PETSC_REAL or PETSC_SCALAR.
323: If the field is complex and datatype is PETSC_REAL, the real part of the
324: field is returned.
327: Output Parameter:
328: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
329: If B is not NULL, the values of the field are written in this array, varying first by component,
330: then by point.
331: . D - pointer to data of size d * c * n * sizeof(datatype).
332: If D is not NULL, the values of the field's spatial derivatives are written in this array,
333: varying first by the partial derivative component, then by field component, then by point.
334: - H - pointer to data of size d * d * c * n * sizeof(datatype).
335: If H is not NULL, the values of the field's second spatial derivatives are written in this array,
336: varying first by the second partial derivative component, then by field component, then by point.
338: Level: intermediate
340: .seealso: DMFieldGetNumComponents(), DMFieldEvaluate(), DMFieldEvaluateFE()
341: @*/
342: PetscErrorCode DMFieldEvaluateFV(DMField field, IS cellIS, PetscDataType datatype, void *B, void *D, void *H)
343: {
352: if (field->ops->evaluateFV) {
353: (*field->ops->evaluateFV) (field, cellIS, datatype, B, D, H);
354: } else SETERRQ (PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented for this type");
355: return(0);
356: }
358: /*@
359: DMFieldGetDegree - Get the polynomial degree of a field when pulled back onto the
360: reference element
362: Not collective
364: Input Arguments:
365: + field - the DMField object
366: - cellIS - the index set of points over which we want know the invariance
368: Output Arguments:
369: + minDegree - the degree of the largest polynomial space contained in the field on each element
370: - maxDegree - the largest degree of the smallest polynomial space containing the field on any element
372: Level: intermediate
374: .seealso: DMFieldEvaluateFE()
375: @*/
376: PetscErrorCode DMFieldGetDegree(DMField field, IS cellIS, PetscInt *minDegree, PetscInt *maxDegree)
377: {
386: if (minDegree) *minDegree = -1;
387: if (maxDegree) *maxDegree = PETSC_MAX_INT;
389: if (field->ops->getDegree) {
390: (*field->ops->getDegree) (field,cellIS,minDegree,maxDegree);
391: }
392: return(0);
393: }
395: /*@
396: DMFieldCreateDefaultQuadrature - Creates a quadrature sufficient to integrate the field on the selected
397: points via pullback onto the reference element
399: Not collective
401: Input Arguments:
402: + field - the DMField object
403: - pointIS - the index set of points over which we wish to integrate the field
405: Output Arguments:
406: . quad - a PetscQuadrature object
408: Level: developer
410: .seealso: DMFieldEvaluteFE(), DMFieldGetDegree()
411: @*/
412: PetscErrorCode DMFieldCreateDefaultQuadrature(DMField field, IS pointIS, PetscQuadrature *quad)
413: {
421: *quad = NULL;
422: if (field->ops->createDefaultQuadrature) {
423: (*field->ops->createDefaultQuadrature)(field, pointIS, quad);
424: }
425: return(0);
426: }
428: /*@C
429: DMFieldCreateFEGeom - Compute and create the geometric factors of a coordinate field
431: Not collective
433: Input Arguments:
434: + field - the DMField object
435: . pointIS - the index set of points over which we wish to integrate the field
436: . quad - the quadrature points at which to evaluate the geometric factors
437: - faceData - whether additional data for facets (the normal vectors and adjacent cells) should
438: be calculated
440: Output Arguments:
441: . geom - the geometric factors
443: Level: developer
445: .seealso: DMFieldEvaluateFE(), DMFieldCreateDefaulteQuadrature(), DMFieldGetDegree()
446: @*/
447: PetscErrorCode DMFieldCreateFEGeom(DMField field, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
448: {
449: PetscInt dim, dE;
450: PetscInt nPoints;
451: PetscInt maxDegree;
452: PetscFEGeom *g;
459: ISGetLocalSize(pointIS,&nPoints);
460: dE = field->numComponents;
461: PetscFEGeomCreate(quad,nPoints,dE,faceData,&g);
462: DMFieldEvaluateFE(field,pointIS,quad,PETSC_REAL,g->v,g->J,NULL);
463: dim = g->dim;
464: if (dE > dim) {
465: /* space out J and make square Jacobians */
466: PetscInt i, j, k, N = g->numPoints * g->numCells;
468: for (i = N-1; i >= 0; i--) {
469: PetscReal J[9];
471: for (j = 0; j < dE; j++) {
472: for (k = 0; k < dim; k++) {
473: J[j*dE + k] = g->J[i*dE*dim + j*dim + k];
474: }
475: }
476: switch (dim) {
477: case 0:
478: for (j = 0; j < dE; j++) {
479: for (k = 0; k < dE; k++) {
480: J[j * dE + k] = (j == k) ? 1. : 0.;
481: }
482: }
483: break;
484: case 1:
485: if (dE == 2) {
486: PetscReal norm = PetscSqrtReal(J[0] * J[0] + J[2] * J[2]);
488: J[1] = -J[2] / norm;
489: J[3] = J[0] / norm;
490: } else {
491: PetscReal inorm = 1./PetscSqrtReal(J[0] * J[0] + J[3] * J[3] + J[6] * J[6]);
492: PetscReal x = J[0] * inorm;
493: PetscReal y = J[3] * inorm;
494: PetscReal z = J[6] * inorm;
496: if (x > 0.) {
497: PetscReal inv1pX = 1./ (1. + x);
499: J[1] = -y; J[2] = -z;
500: J[4] = 1. - y*y*inv1pX; J[5] = -y*z*inv1pX;
501: J[7] = -y*z*inv1pX; J[8] = 1. - z*z*inv1pX;
502: } else {
503: PetscReal inv1mX = 1./ (1. - x);
505: J[1] = z; J[2] = y;
506: J[4] = -y*z*inv1mX; J[5] = 1. - y*y*inv1mX;
507: J[7] = 1. - z*z*inv1mX; J[8] = -y*z*inv1mX;
508: }
509: }
510: break;
511: case 2:
512: {
513: PetscReal inorm;
515: J[2] = J[3] * J[7] - J[6] * J[4];
516: J[5] = J[6] * J[1] - J[0] * J[7];
517: J[8] = J[0] * J[4] - J[3] * J[1];
519: inorm = 1./ PetscSqrtReal(J[2]*J[2] + J[5]*J[5] + J[8]*J[8]);
521: J[2] *= inorm;
522: J[5] *= inorm;
523: J[8] *= inorm;
524: }
525: break;
526: }
527: for (j = 0; j < dE*dE; j++) {
528: g->J[i*dE*dE + j] = J[j];
529: }
530: }
531: }
532: PetscFEGeomComplete(g);
533: DMFieldGetDegree(field,pointIS,NULL,&maxDegree);
534: g->isAffine = (maxDegree <= 1) ? PETSC_TRUE : PETSC_FALSE;
535: if (faceData) {
536: if (!field->ops->computeFaceData) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "DMField implementation does not compute face data\n");
537: (*field->ops->computeFaceData) (field, pointIS, quad, g);
538: }
539: *geom = g;
540: return(0);
541: }