Actual source code: dmfieldshell.c
petsc-3.14.6 2021-03-30
1: #include <petsc/private/dmfieldimpl.h>
3: typedef struct _n_DMField_Shell
4: {
5: void *ctx;
6: PetscErrorCode (*destroy) (DMField);
7: }
8: DMField_Shell;
10: PetscErrorCode DMFieldShellGetContext(DMField field, void *ctx)
11: {
12: PetscBool flg;
18: PetscObjectTypeCompare((PetscObject)field,DMFIELDSHELL,&flg);
19: if (flg) *(void**)ctx = ((DMField_Shell*)(field->data))->ctx;
20: else SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Cannot get context from non-shell shield");
21: return(0);
22: }
24: static PetscErrorCode DMFieldDestroy_Shell(DMField field)
25: {
26: DMField_Shell *shell = (DMField_Shell *) field->data;
30: if (shell->destroy) {(*(shell->destroy)) (field);}
31: PetscFree(field->data);
32: return(0);
33: }
35: PetscErrorCode DMFieldShellEvaluateFEDefault(DMField field, IS pointIS, PetscQuadrature quad, PetscDataType type, void *B, void *D, void *H)
36: {
37: DM dm = field->dm;
38: DMField coordField;
39: PetscFEGeom *geom;
40: Vec pushforward;
41: PetscInt dimC, dim, numPoints, Nq, p, Nc;
42: PetscScalar *pfArray;
43: PetscErrorCode ierr;
46: Nc = field->numComponents;
47: DMGetCoordinateField(dm, &coordField);
48: DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FALSE, &geom);
49: DMGetCoordinateDim(dm, &dimC);
50: PetscQuadratureGetData(quad, &dim, NULL, &Nq, NULL, NULL);
51: ISGetLocalSize(pointIS, &numPoints);
52: PetscMalloc1(dimC * Nq * numPoints, &pfArray);
53: for (p = 0; p < numPoints * dimC * Nq; p++) pfArray[p] = (PetscScalar) geom->v[p];
54: VecCreateMPIWithArray(PetscObjectComm((PetscObject)pointIS), dimC, dimC * Nq * numPoints, PETSC_DETERMINE, pfArray, &pushforward);
55: DMFieldEvaluate(field, pushforward, type, B, D, H);
56: /* TODO: handle covariant/contravariant pullbacks */
57: if (D) {
58: if (type == PETSC_SCALAR) {
59: PetscScalar *sD = (PetscScalar *) D;
60: PetscInt q;
62: for (p = 0; p < numPoints * Nq; p++) {
63: for (q = 0; q < Nc; q++) {
64: PetscScalar d[3];
66: PetscInt i, j;
68: for (i = 0; i < dimC; i++) d[i] = sD[(p * Nc + q) * dimC + i];
69: for (i = 0; i < dimC; i++) sD[(p * Nc + q) * dimC + i] = 0.;
70: for (i = 0; i < dimC; i++) {
71: for (j = 0; j < dimC; j++) {
72: sD[(p * Nc + q) * dimC + i] += geom->J[(p * dimC + j) * dimC + i] * d[j];
73: }
74: }
75: }
76: }
77: } else {
78: PetscReal *rD = (PetscReal *) D;
79: PetscInt q;
81: for (p = 0; p < numPoints * Nq; p++) {
82: for (q = 0; q < Nc; q++) {
83: PetscReal d[3];
85: PetscInt i, j;
87: for (i = 0; i < dimC; i++) d[i] = rD[(p * Nc + q) * dimC + i];
88: for (i = 0; i < dimC; i++) rD[(p * Nc + q) * dimC + i] = 0.;
89: for (i = 0; i < dimC; i++) {
90: for (j = 0; j < dimC; j++) {
91: rD[(p * Nc + q) * dimC + i] += geom->J[(p * dimC + j) * dimC + i] * d[j];
92: }
93: }
94: }
95: }
96: }
97: }
98: if (H) {
99: if (type == PETSC_SCALAR) {
100: PetscScalar *sH = (PetscScalar *) H;
101: PetscInt q;
103: for (p = 0; p < numPoints * Nq; p++) {
104: for (q = 0; q < Nc; q++) {
105: PetscScalar d[3][3];
107: PetscInt i, j, k, l;
109: for (i = 0; i < dimC; i++) for (j = 0; j < dimC; j++) d[i][j] = sH[((p * Nc + q) * dimC + i) * dimC + j];
110: for (i = 0; i < dimC; i++) for (j = 0; j < dimC; j++) sH[((p * Nc + q) * dimC + i) * dimC + j] = 0.;
111: for (i = 0; i < dimC; i++) {
112: for (j = 0; j < dimC; j++) {
113: for (k = 0; k < dimC; k++) {
114: for (l = 0; l < dimC; l++) {
115: sH[((p * Nc + q) * dimC + i) * dimC + j] += geom->J[(p * dimC + k) * dimC + i] * geom->J[(p * dimC + l) * dimC + j] * d[k][l];
116: }
117: }
118: }
119: }
120: }
121: }
122: } else {
123: PetscReal *rH = (PetscReal *) H;
124: PetscInt q;
126: for (p = 0; p < numPoints * Nq; p++) {
127: for (q = 0; q < Nc; q++) {
128: PetscReal d[3][3];
130: PetscInt i, j, k, l;
132: for (i = 0; i < dimC; i++) for (j = 0; j < dimC; j++) d[i][j] = rH[((p * Nc + q) * dimC + i) * dimC + j];
133: for (i = 0; i < dimC; i++) for (j = 0; j < dimC; j++) rH[((p * Nc + q) * dimC + i) * dimC + j] = 0.;
134: for (i = 0; i < dimC; i++) {
135: for (j = 0; j < dimC; j++) {
136: for (k = 0; k < dimC; k++) {
137: for (l = 0; l < dimC; l++) {
138: rH[((p * Nc + q) * dimC + i) * dimC + j] += geom->J[(p * dimC + k) * dimC + i] * geom->J[(p * dimC + l) * dimC + j] * d[k][l];
139: }
140: }
141: }
142: }
143: }
144: }
145: }
146: }
147: VecDestroy(&pushforward);
148: PetscFree(pfArray);
149: PetscFEGeomDestroy(&geom);
150: return(0);
151: }
153: PetscErrorCode DMFieldShellEvaluateFVDefault(DMField field, IS pointIS, PetscDataType type, void *B, void *D, void *H)
154: {
155: DM dm = field->dm;
156: DMField coordField;
157: PetscFEGeom *geom;
158: Vec pushforward;
159: PetscInt dimC, dim, numPoints, Nq, p;
160: PetscScalar *pfArray;
161: PetscQuadrature quad;
162: PetscErrorCode ierr;
165: DMGetCoordinateField(dm, &coordField);
166: DMFieldCreateDefaultQuadrature(coordField, pointIS, &quad);
167: if (!quad) SETERRQ(PetscObjectComm((PetscObject) pointIS), PETSC_ERR_ARG_WRONGSTATE, "coordinate field must have default quadrature for FV computation");
168: DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FALSE, &geom);
169: DMGetCoordinateDim(dm, &dimC);
170: PetscQuadratureGetData(quad, &dim, NULL, &Nq, NULL, NULL);
171: if (Nq != 1) SETERRQ(PetscObjectComm((PetscObject) quad), PETSC_ERR_ARG_WRONGSTATE, "quadrature must have only one point");
172: ISGetLocalSize(pointIS, &numPoints);
173: PetscMalloc1(dimC * numPoints, &pfArray);
174: for (p = 0; p < numPoints * dimC; p++) pfArray[p] = (PetscScalar) geom->v[p];
175: VecCreateMPIWithArray(PetscObjectComm((PetscObject)pointIS), dimC, dimC * numPoints, PETSC_DETERMINE, pfArray, &pushforward);
176: DMFieldEvaluate(field, pushforward, type, B, D, H);
177: PetscQuadratureDestroy(&quad);
178: VecDestroy(&pushforward);
179: PetscFree(pfArray);
180: PetscFEGeomDestroy(&geom);
181: return(0);
182: }
184: PetscErrorCode DMFieldShellSetDestroy(DMField field, PetscErrorCode (*destroy)(DMField))
185: {
186: DMField_Shell *shell = (DMField_Shell *) field->data;
190: shell->destroy = destroy;
191: return(0);
192: }
194: PetscErrorCode DMFieldShellSetEvaluate(DMField field, PetscErrorCode (*evaluate)(DMField,Vec,PetscDataType,void*,void*,void*))
195: {
198: field->ops->evaluate = evaluate;
199: return(0);
200: }
202: PetscErrorCode DMFieldShellSetEvaluateFE(DMField field, PetscErrorCode (*evaluateFE)(DMField,IS,PetscQuadrature,PetscDataType,void*,void*,void*))
203: {
206: field->ops->evaluateFE = evaluateFE;
207: return(0);
208: }
210: PetscErrorCode DMFieldShellSetEvaluateFV(DMField field, PetscErrorCode (*evaluateFV)(DMField,IS,PetscDataType,void*,void*,void*))
211: {
214: field->ops->evaluateFV = evaluateFV;
215: return(0);
216: }
218: PetscErrorCode DMFieldShellSetGetDegree(DMField field, PetscErrorCode (*getDegree)(DMField,IS,PetscInt*,PetscInt*))
219: {
222: field->ops->getDegree = getDegree;
223: return(0);
224: }
226: PetscErrorCode DMFieldShellSetCreateDefaultQuadrature(DMField field, PetscErrorCode (*createDefaultQuadrature)(DMField,IS,PetscQuadrature*))
227: {
230: field->ops->createDefaultQuadrature = createDefaultQuadrature;
231: return(0);
232: }
234: static PetscErrorCode DMFieldInitialize_Shell(DMField field)
235: {
237: field->ops->destroy = DMFieldDestroy_Shell;
238: field->ops->evaluate = NULL;
239: field->ops->evaluateFE = DMFieldShellEvaluateFEDefault;
240: field->ops->evaluateFV = DMFieldShellEvaluateFVDefault;
241: field->ops->getDegree = NULL;
242: field->ops->createDefaultQuadrature = NULL;
243: field->ops->view = NULL;
244: return(0);
245: }
247: PETSC_INTERN PetscErrorCode DMFieldCreate_Shell(DMField field)
248: {
249: DMField_Shell *shell;
253: PetscNewLog(field,&shell);
254: field->data = shell;
255: DMFieldInitialize_Shell(field);
256: return(0);
257: }
259: PetscErrorCode DMFieldCreateShell(DM dm, PetscInt numComponents, DMFieldContinuity continuity, void *ctx, DMField *field)
260: {
261: DMField b;
262: DMField_Shell *shell;
269: DMFieldCreate(dm,numComponents,continuity,&b);
270: DMFieldSetType(b,DMFIELDSHELL);
271: shell = (DMField_Shell *) b->data;
272: shell->ctx = ctx;
273: *field = b;
274: return(0);
275: }