Actual source code: dmfieldshell.c
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;
16: PetscObjectTypeCompare((PetscObject)field,DMFIELDSHELL,&flg);
17: if (flg) *(void**)ctx = ((DMField_Shell*)(field->data))->ctx;
18: else SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Cannot get context from non-shell shield");
19: return 0;
20: }
22: static PetscErrorCode DMFieldDestroy_Shell(DMField field)
23: {
24: DMField_Shell *shell = (DMField_Shell *) field->data;
26: if (shell->destroy) (*(shell->destroy)) (field);
27: PetscFree(field->data);
28: return 0;
29: }
31: PetscErrorCode DMFieldShellEvaluateFEDefault(DMField field, IS pointIS, PetscQuadrature quad, PetscDataType type, void *B, void *D, void *H)
32: {
33: DM dm = field->dm;
34: DMField coordField;
35: PetscFEGeom *geom;
36: Vec pushforward;
37: PetscInt dimC, dim, numPoints, Nq, p, Nc;
38: PetscScalar *pfArray;
40: Nc = field->numComponents;
41: DMGetCoordinateField(dm, &coordField);
42: DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FALSE, &geom);
43: DMGetCoordinateDim(dm, &dimC);
44: PetscQuadratureGetData(quad, &dim, NULL, &Nq, NULL, NULL);
45: ISGetLocalSize(pointIS, &numPoints);
46: PetscMalloc1(dimC * Nq * numPoints, &pfArray);
47: for (p = 0; p < numPoints * dimC * Nq; p++) pfArray[p] = (PetscScalar) geom->v[p];
48: VecCreateMPIWithArray(PetscObjectComm((PetscObject)pointIS), dimC, dimC * Nq * numPoints, PETSC_DETERMINE, pfArray, &pushforward);
49: DMFieldEvaluate(field, pushforward, type, B, D, H);
50: /* TODO: handle covariant/contravariant pullbacks */
51: if (D) {
52: if (type == PETSC_SCALAR) {
53: PetscScalar *sD = (PetscScalar *) D;
54: PetscInt q;
56: for (p = 0; p < numPoints * Nq; p++) {
57: for (q = 0; q < Nc; q++) {
58: PetscScalar d[3];
60: PetscInt i, j;
62: for (i = 0; i < dimC; i++) d[i] = sD[(p * Nc + q) * dimC + i];
63: for (i = 0; i < dimC; i++) sD[(p * Nc + q) * dimC + i] = 0.;
64: for (i = 0; i < dimC; i++) {
65: for (j = 0; j < dimC; j++) {
66: sD[(p * Nc + q) * dimC + i] += geom->J[(p * dimC + j) * dimC + i] * d[j];
67: }
68: }
69: }
70: }
71: } else {
72: PetscReal *rD = (PetscReal *) D;
73: PetscInt q;
75: for (p = 0; p < numPoints * Nq; p++) {
76: for (q = 0; q < Nc; q++) {
77: PetscReal d[3];
79: PetscInt i, j;
81: for (i = 0; i < dimC; i++) d[i] = rD[(p * Nc + q) * dimC + i];
82: for (i = 0; i < dimC; i++) rD[(p * Nc + q) * dimC + i] = 0.;
83: for (i = 0; i < dimC; i++) {
84: for (j = 0; j < dimC; j++) {
85: rD[(p * Nc + q) * dimC + i] += geom->J[(p * dimC + j) * dimC + i] * d[j];
86: }
87: }
88: }
89: }
90: }
91: }
92: if (H) {
93: if (type == PETSC_SCALAR) {
94: PetscScalar *sH = (PetscScalar *) H;
95: PetscInt q;
97: for (p = 0; p < numPoints * Nq; p++) {
98: for (q = 0; q < Nc; q++) {
99: PetscScalar d[3][3];
101: PetscInt i, j, k, l;
103: for (i = 0; i < dimC; i++) for (j = 0; j < dimC; j++) d[i][j] = sH[((p * Nc + q) * dimC + i) * dimC + j];
104: for (i = 0; i < dimC; i++) for (j = 0; j < dimC; j++) sH[((p * Nc + q) * dimC + i) * dimC + j] = 0.;
105: for (i = 0; i < dimC; i++) {
106: for (j = 0; j < dimC; j++) {
107: for (k = 0; k < dimC; k++) {
108: for (l = 0; l < dimC; l++) {
109: 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];
110: }
111: }
112: }
113: }
114: }
115: }
116: } else {
117: PetscReal *rH = (PetscReal *) H;
118: PetscInt q;
120: for (p = 0; p < numPoints * Nq; p++) {
121: for (q = 0; q < Nc; q++) {
122: PetscReal d[3][3];
124: PetscInt i, j, k, l;
126: for (i = 0; i < dimC; i++) for (j = 0; j < dimC; j++) d[i][j] = rH[((p * Nc + q) * dimC + i) * dimC + j];
127: for (i = 0; i < dimC; i++) for (j = 0; j < dimC; j++) rH[((p * Nc + q) * dimC + i) * dimC + j] = 0.;
128: for (i = 0; i < dimC; i++) {
129: for (j = 0; j < dimC; j++) {
130: for (k = 0; k < dimC; k++) {
131: for (l = 0; l < dimC; l++) {
132: 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];
133: }
134: }
135: }
136: }
137: }
138: }
139: }
140: }
141: VecDestroy(&pushforward);
142: PetscFree(pfArray);
143: PetscFEGeomDestroy(&geom);
144: return 0;
145: }
147: PetscErrorCode DMFieldShellEvaluateFVDefault(DMField field, IS pointIS, PetscDataType type, void *B, void *D, void *H)
148: {
149: DM dm = field->dm;
150: DMField coordField;
151: PetscFEGeom *geom;
152: Vec pushforward;
153: PetscInt dimC, dim, numPoints, Nq, p;
154: PetscScalar *pfArray;
155: PetscQuadrature quad;
156: MPI_Comm comm;
158: PetscObjectGetComm((PetscObject) field, &comm);
159: DMGetDimension(dm, &dim);
160: DMGetCoordinateDim(dm, &dimC);
161: DMGetCoordinateField(dm, &coordField);
162: DMFieldGetFVQuadrature_Internal(coordField, pointIS, &quad);
164: PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
166: DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FALSE, &geom);
167: ISGetLocalSize(pointIS, &numPoints);
168: PetscMalloc1(dimC * numPoints, &pfArray);
169: for (p = 0; p < numPoints * dimC; p++) pfArray[p] = (PetscScalar) geom->v[p];
170: VecCreateMPIWithArray(PetscObjectComm((PetscObject)pointIS), dimC, dimC * numPoints, PETSC_DETERMINE, pfArray, &pushforward);
171: DMFieldEvaluate(field, pushforward, type, B, D, H);
172: PetscQuadratureDestroy(&quad);
173: VecDestroy(&pushforward);
174: PetscFree(pfArray);
175: PetscFEGeomDestroy(&geom);
176: return 0;
177: }
179: PetscErrorCode DMFieldShellSetDestroy(DMField field, PetscErrorCode (*destroy)(DMField))
180: {
181: DMField_Shell *shell = (DMField_Shell *) field->data;
184: shell->destroy = destroy;
185: return 0;
186: }
188: PetscErrorCode DMFieldShellSetEvaluate(DMField field, PetscErrorCode (*evaluate)(DMField,Vec,PetscDataType,void*,void*,void*))
189: {
191: field->ops->evaluate = evaluate;
192: return 0;
193: }
195: PetscErrorCode DMFieldShellSetEvaluateFE(DMField field, PetscErrorCode (*evaluateFE)(DMField,IS,PetscQuadrature,PetscDataType,void*,void*,void*))
196: {
198: field->ops->evaluateFE = evaluateFE;
199: return 0;
200: }
202: PetscErrorCode DMFieldShellSetEvaluateFV(DMField field, PetscErrorCode (*evaluateFV)(DMField,IS,PetscDataType,void*,void*,void*))
203: {
205: field->ops->evaluateFV = evaluateFV;
206: return 0;
207: }
209: PetscErrorCode DMFieldShellSetGetDegree(DMField field, PetscErrorCode (*getDegree)(DMField,IS,PetscInt*,PetscInt*))
210: {
212: field->ops->getDegree = getDegree;
213: return 0;
214: }
216: PetscErrorCode DMFieldShellSetCreateDefaultQuadrature(DMField field, PetscErrorCode (*createDefaultQuadrature)(DMField,IS,PetscQuadrature*))
217: {
219: field->ops->createDefaultQuadrature = createDefaultQuadrature;
220: return 0;
221: }
223: static PetscErrorCode DMFieldInitialize_Shell(DMField field)
224: {
225: field->ops->destroy = DMFieldDestroy_Shell;
226: field->ops->evaluate = NULL;
227: field->ops->evaluateFE = DMFieldShellEvaluateFEDefault;
228: field->ops->evaluateFV = DMFieldShellEvaluateFVDefault;
229: field->ops->getDegree = NULL;
230: field->ops->createDefaultQuadrature = NULL;
231: field->ops->view = NULL;
232: return 0;
233: }
235: PETSC_INTERN PetscErrorCode DMFieldCreate_Shell(DMField field)
236: {
237: DMField_Shell *shell;
239: PetscNewLog(field,&shell);
240: field->data = shell;
241: DMFieldInitialize_Shell(field);
242: return 0;
243: }
245: PetscErrorCode DMFieldCreateShell(DM dm, PetscInt numComponents, DMFieldContinuity continuity, void *ctx, DMField *field)
246: {
247: DMField b;
248: DMField_Shell *shell;
253: DMFieldCreate(dm,numComponents,continuity,&b);
254: DMFieldSetType(b,DMFIELDSHELL);
255: shell = (DMField_Shell *) b->data;
256: shell->ctx = ctx;
257: *field = b;
258: return 0;
259: }