Actual source code: dmfieldshell.c
1: #include <petsc/private/dmfieldimpl.h>
3: typedef struct _n_DMField_Shell {
4: void *ctx;
5: PetscErrorCode (*destroy)(DMField);
6: } DMField_Shell;
8: PetscErrorCode DMFieldShellGetContext(DMField field, void *ctx)
9: {
10: PetscBool flg;
14: PetscObjectTypeCompare((PetscObject)field, DMFIELDSHELL, &flg);
15: if (flg) *(void **)ctx = ((DMField_Shell *)(field->data))->ctx;
16: else SETERRQ(PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "Cannot get context from non-shell shield");
17: return 0;
18: }
20: static PetscErrorCode DMFieldDestroy_Shell(DMField field)
21: {
22: DMField_Shell *shell = (DMField_Shell *)field->data;
24: if (shell->destroy) (*(shell->destroy))(field);
25: PetscFree(field->data);
26: return 0;
27: }
29: PetscErrorCode DMFieldShellEvaluateFEDefault(DMField field, IS pointIS, PetscQuadrature quad, PetscDataType type, void *B, void *D, void *H)
30: {
31: DM dm = field->dm;
32: DMField coordField;
33: PetscFEGeom *geom;
34: Vec pushforward;
35: PetscInt dimC, dim, numPoints, Nq, p, Nc;
36: PetscScalar *pfArray;
38: Nc = field->numComponents;
39: DMGetCoordinateField(dm, &coordField);
40: DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FALSE, &geom);
41: DMGetCoordinateDim(dm, &dimC);
42: PetscQuadratureGetData(quad, &dim, NULL, &Nq, NULL, NULL);
43: ISGetLocalSize(pointIS, &numPoints);
44: PetscMalloc1(dimC * Nq * numPoints, &pfArray);
45: for (p = 0; p < numPoints * dimC * Nq; p++) pfArray[p] = (PetscScalar)geom->v[p];
46: VecCreateMPIWithArray(PetscObjectComm((PetscObject)pointIS), dimC, dimC * Nq * numPoints, PETSC_DETERMINE, pfArray, &pushforward);
47: DMFieldEvaluate(field, pushforward, type, B, D, H);
48: /* TODO: handle covariant/contravariant pullbacks */
49: if (D) {
50: if (type == PETSC_SCALAR) {
51: PetscScalar *sD = (PetscScalar *)D;
52: PetscInt q;
54: for (p = 0; p < numPoints * Nq; p++) {
55: for (q = 0; q < Nc; q++) {
56: PetscScalar d[3];
58: PetscInt i, j;
60: for (i = 0; i < dimC; i++) d[i] = sD[(p * Nc + q) * dimC + i];
61: for (i = 0; i < dimC; i++) sD[(p * Nc + q) * dimC + i] = 0.;
62: for (i = 0; i < dimC; i++) {
63: for (j = 0; j < dimC; j++) sD[(p * Nc + q) * dimC + i] += geom->J[(p * dimC + j) * dimC + i] * d[j];
64: }
65: }
66: }
67: } else {
68: PetscReal *rD = (PetscReal *)D;
69: PetscInt q;
71: for (p = 0; p < numPoints * Nq; p++) {
72: for (q = 0; q < Nc; q++) {
73: PetscReal d[3];
75: PetscInt i, j;
77: for (i = 0; i < dimC; i++) d[i] = rD[(p * Nc + q) * dimC + i];
78: for (i = 0; i < dimC; i++) rD[(p * Nc + q) * dimC + i] = 0.;
79: for (i = 0; i < dimC; i++) {
80: for (j = 0; j < dimC; j++) rD[(p * Nc + q) * dimC + i] += geom->J[(p * dimC + j) * dimC + i] * d[j];
81: }
82: }
83: }
84: }
85: }
86: if (H) {
87: if (type == PETSC_SCALAR) {
88: PetscScalar *sH = (PetscScalar *)H;
89: PetscInt q;
91: for (p = 0; p < numPoints * Nq; p++) {
92: for (q = 0; q < Nc; q++) {
93: PetscScalar d[3][3];
95: PetscInt i, j, k, l;
97: for (i = 0; i < dimC; i++)
98: for (j = 0; j < dimC; j++) d[i][j] = sH[((p * Nc + q) * dimC + i) * dimC + j];
99: for (i = 0; i < dimC; i++)
100: for (j = 0; j < dimC; j++) sH[((p * Nc + q) * dimC + i) * dimC + j] = 0.;
101: for (i = 0; i < dimC; i++) {
102: for (j = 0; j < dimC; j++) {
103: for (k = 0; k < dimC; k++) {
104: for (l = 0; l < dimC; l++) 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];
105: }
106: }
107: }
108: }
109: }
110: } else {
111: PetscReal *rH = (PetscReal *)H;
112: PetscInt q;
114: for (p = 0; p < numPoints * Nq; p++) {
115: for (q = 0; q < Nc; q++) {
116: PetscReal d[3][3];
118: PetscInt i, j, k, l;
120: for (i = 0; i < dimC; i++)
121: for (j = 0; j < dimC; j++) d[i][j] = rH[((p * Nc + q) * dimC + i) * dimC + j];
122: for (i = 0; i < dimC; i++)
123: for (j = 0; j < dimC; j++) rH[((p * Nc + q) * dimC + i) * dimC + j] = 0.;
124: for (i = 0; i < dimC; i++) {
125: for (j = 0; j < dimC; j++) {
126: for (k = 0; k < dimC; k++) {
127: for (l = 0; l < dimC; l++) 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];
128: }
129: }
130: }
131: }
132: }
133: }
134: }
135: VecDestroy(&pushforward);
136: PetscFree(pfArray);
137: PetscFEGeomDestroy(&geom);
138: return 0;
139: }
141: PetscErrorCode DMFieldShellEvaluateFVDefault(DMField field, IS pointIS, PetscDataType type, void *B, void *D, void *H)
142: {
143: DM dm = field->dm;
144: DMField coordField;
145: PetscFEGeom *geom;
146: Vec pushforward;
147: PetscInt dimC, dim, numPoints, Nq, p;
148: PetscScalar *pfArray;
149: PetscQuadrature quad;
150: MPI_Comm comm;
152: PetscObjectGetComm((PetscObject)field, &comm);
153: DMGetDimension(dm, &dim);
154: DMGetCoordinateDim(dm, &dimC);
155: DMGetCoordinateField(dm, &coordField);
156: DMFieldGetFVQuadrature_Internal(coordField, pointIS, &quad);
158: PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
160: DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FALSE, &geom);
161: ISGetLocalSize(pointIS, &numPoints);
162: PetscMalloc1(dimC * numPoints, &pfArray);
163: for (p = 0; p < numPoints * dimC; p++) pfArray[p] = (PetscScalar)geom->v[p];
164: VecCreateMPIWithArray(PetscObjectComm((PetscObject)pointIS), dimC, dimC * numPoints, PETSC_DETERMINE, pfArray, &pushforward);
165: DMFieldEvaluate(field, pushforward, type, B, D, H);
166: PetscQuadratureDestroy(&quad);
167: VecDestroy(&pushforward);
168: PetscFree(pfArray);
169: PetscFEGeomDestroy(&geom);
170: return 0;
171: }
173: PetscErrorCode DMFieldShellSetDestroy(DMField field, PetscErrorCode (*destroy)(DMField))
174: {
175: DMField_Shell *shell = (DMField_Shell *)field->data;
178: shell->destroy = destroy;
179: return 0;
180: }
182: PetscErrorCode DMFieldShellSetEvaluate(DMField field, PetscErrorCode (*evaluate)(DMField, Vec, PetscDataType, void *, void *, void *))
183: {
185: field->ops->evaluate = evaluate;
186: return 0;
187: }
189: PetscErrorCode DMFieldShellSetEvaluateFE(DMField field, PetscErrorCode (*evaluateFE)(DMField, IS, PetscQuadrature, PetscDataType, void *, void *, void *))
190: {
192: field->ops->evaluateFE = evaluateFE;
193: return 0;
194: }
196: PetscErrorCode DMFieldShellSetEvaluateFV(DMField field, PetscErrorCode (*evaluateFV)(DMField, IS, PetscDataType, void *, void *, void *))
197: {
199: field->ops->evaluateFV = evaluateFV;
200: return 0;
201: }
203: PetscErrorCode DMFieldShellSetGetDegree(DMField field, PetscErrorCode (*getDegree)(DMField, IS, PetscInt *, PetscInt *))
204: {
206: field->ops->getDegree = getDegree;
207: return 0;
208: }
210: PetscErrorCode DMFieldShellSetCreateDefaultQuadrature(DMField field, PetscErrorCode (*createDefaultQuadrature)(DMField, IS, PetscQuadrature *))
211: {
213: field->ops->createDefaultQuadrature = createDefaultQuadrature;
214: return 0;
215: }
217: static PetscErrorCode DMFieldInitialize_Shell(DMField field)
218: {
219: field->ops->destroy = DMFieldDestroy_Shell;
220: field->ops->evaluate = NULL;
221: field->ops->evaluateFE = DMFieldShellEvaluateFEDefault;
222: field->ops->evaluateFV = DMFieldShellEvaluateFVDefault;
223: field->ops->getDegree = NULL;
224: field->ops->createDefaultQuadrature = NULL;
225: field->ops->view = NULL;
226: return 0;
227: }
229: PETSC_INTERN PetscErrorCode DMFieldCreate_Shell(DMField field)
230: {
231: DMField_Shell *shell;
233: PetscNew(&shell);
234: field->data = shell;
235: DMFieldInitialize_Shell(field);
236: return 0;
237: }
239: PetscErrorCode DMFieldCreateShell(DM dm, PetscInt numComponents, DMFieldContinuity continuity, void *ctx, DMField *field)
240: {
241: DMField b;
242: DMField_Shell *shell;
247: DMFieldCreate(dm, numComponents, continuity, &b);
248: DMFieldSetType(b, DMFIELDSHELL);
249: shell = (DMField_Shell *)b->data;
250: shell->ctx = ctx;
251: *field = b;
252: return 0;
253: }